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Prediction  of  Physical  Analyte  Data  Based  on  Interactions  with  a  Chemical  Vapor 
Detection  System 

Introduction 

The  use  of  Quantitative  Structure-Activity  Relationship  (QSAR)  and  Quantitative  Structure-Property 
Relationship  (QSPR)  techniques  has  allowed  for  the  characterization  of  untested  molecules  based  on 
similarities  to  a  library  of  experimentally  characterized  molecules.  Essentially,  the  QSAR/QSPR 
method  involves  the  development  of  a  prediction  model  that  matches  an  array  of  experimental  data. 
This  model  can  be  of  effectively  any  quantitative  nature,  and  many  methods  have  been  used,  such  as 
multiple  linear  regression,  partial  least  squares,  neural  network  analysis,  and  genetic  functional 
analysis. 

QSAR/QSPR  has  been  particularly  useful  for  drug  discovery,  for  which  a  fast  method  of  pre¬ 
screening  drug  candidates  based  on  expected  characteristics  is  necessary.  Additionally,  this 
technique  has  been  used  to  predict  fundamental  properties  of  molecules,  such  as  thermodynamic 
constants,  based  on  their  structures.  However,  these  methods  often  employ  semi-quantitative 
"features,"  such  as  the  number  of  double  bonds  or  benzene  ring  in  a  structure1,  or  other 
thermodynamic  data2.  However,  such  data  can  be  subjective  (in  the  case  of  semi-quantitative 
features)  or  not  available  for  all  possible  molecules  of  interest  (in  the  case  of  thermodynamic  data). 
For  these  reasons,  a  method  that  is  based  on  simple,  fast  experimental  methods  may  in  some 
situations  be  desirable. 

Recent  studies  have  shown  that  the  responses  of  polymer-based  chemical  vapor  detectors  can  be 
modeled  by  thermodynamic  relationships  between  the  polymer  sensor  and  analyte  vapor.3  This  is 
reasonable,  considering  that  detector  responses  are  largely  determined  by  the  partition  coefficient  K, 
which  varies  for  every  combination  of  analyte  and  polymer.3,4  Therefore,  it  is  reasonable  that  the 
reverse  should  also  be  possible  -  that  is,  the  responses  of  an  array  of  different  detectors  to  a  library  of 
analytes  should  allow  for  the  development  of  a  predictive  model  that  characterizes  analytes  based 
upon  their  responses  to  a  chemical  detector  array. 

Such  an  attempt  was  made  by  Vaid  et  al.,  in  which  cytochrome-p450  activity  for  a  series  of  alcohols 
was  predicted  based  upon  the  responses  of  these  alcohols  to  a  detector  array.5  In  this  case,  Genetic 
Functional  Analysis  (GFA)6  was  used  to  develop  a  predictive  model  relating  cytochrome-p450 
activity  to  sensor  responses.  GFA  yielded  100  different  models,  which  were  ranked  based  on  their 
lack-of-fit  (LOF,  equation  1). 

eq  1 .  LOF=LSE/(  1  -((c+dp)/M)2) 

Ultimately,  the  equation  used  to  model  cytochrome-p450  activity  for  alcohols  yielded  a  high  degree 
of  statistical  significance  as  determined  by  its  LOF  and  F-test  values.  Because  the  set  was  quite 
small,  cross-validation  was  not  performed  -  rather,  an  estimated  cross-validation  R2  value  was 
determined  from  the  basis  dataset. 


Experimental 


For  our  current  work,  we  seek  to  explore  whether  other  characteristic  analyte  properties  may  be 
predicted  from  sensor  responses.  We  have  explored  bulk  analyte  properties  (water/octanol  partition 
coefficients  (KoW)),  size-related  physical  parameters  (van  der  Waals  area,  volume,  and  radius  of 
gyration),  and  energetic  parameters  such  as  solubility  parameter  ((AH-RT)/Vm)1/2  and  Hansen 
energies,  which  are  decomposed  from  solubility  parameter  as  shown  in  eq  2.: 

Eq  2.  SP  =  ([(HanE)2  +  (HanD)2  +  (HanHB)2]/  Vm)1/2 

Here,  HanE,  HanD,  and  HanHB  are  the  Hansen  electrostatic,  dipole,  and  hydrogen  bonding 
components  of  the  Hildebrandt  cohesive  energy,  respectively;  SP  is  the  solubility  parameter  term, 
and  Vm  is  the  liquid  molar  volume.  Our  research  did  not  focus  on  the  Hansen  hydrogen  bonding 
term  because  of  the  small  number  of  analytes  in  our  set  that  possessed  non-zero  values  for  this 
property.  All  values  were  acquired  from  the  DIPPR  online  database.7  Values  of  KqW  and  solubility 
parameters  were  determined  experimentally,  while  other  values  were  derived  or  calculated. 

An  array  of  40  sensors,  2  copies  each  of  20  different  detectors,  was  exposed  10  times  to  each  of  75 
different  analytes,  which  fall  roughly  into  five  classes:  alcohols,  halides,  aromatics,  hydrocarbons, 
and  esters.  Responses  of  each  sensor  to  each  analyte  exposure  were  determined  by  calculating  the 
fractional  steady-state  resistance  change  of  each  sensor,  ARcq/Rb.  More  details  regarding  this  sensor 
array  and  analyte  set,  used  for  other  work,  are  described  elsewhere.8 

To  form  a  basis  for  prediction  models,  the  data  were  divided  into  a  training  set  (consisting  of  18 
analytes)  which  was  used  to  build  the  models,  a  validation  set  (consisting  of  19  analytes)  used  to 
choose  the  best  model,  and  a  test  set  (consisting  of  37  analytes)  that  was  used  to  evaluate  the  model. 
To  ensure  that  the  test  set  was  interpolating  rather  than  extrapolating  based  on  the  model,  any 
analytes  that  possessed  either  the  highest  or  lowest  values  for  any  property  considered  were  assigned 
to  the  training  set.  The  other  analytes  were  randomly  assigned  set  membership,  while  assuring  that 
each  analyte  class  contributed  roughly  evenly  to  ensure  that  each  class  was  fairly  represented  in  each 
set. 

Previous  work,  both  with  cytochrome-p450  activity  prediction  and  with  prediction  of  the  properties 
here,  was  performed  using  GFA  models  that  utilized  a  handful  of  sensors.  However,  this  method  has 
the  disadvantage  that  a  great  deal  of  useful  information  may  be  contained  in  the  sensors  rejected, 
while  using  more  sensors  would  slow  the  algorithm  and  almost  certainly  result  in  overtraining.  For 
this  reason,  rather  than  choosing  sensors  the  data  was  first  transformed  using  Fisher’s  Linear 
Discriminant  (FLD).9  FLD  takes  as  input  our  40-dimensional  sensor  data,  and  transforms  the  data 
such  that  the  greatest  portion  of  the  data  which  distinguishes  the  different  analytes  is  contained  in  the 
first  few  dimensions.  In  this  way,  well  over  99%  of  the  40-dimensional  information  can  be 
contained  in  as  few  as  five  dimensions.  This  simplifies  the  models  and  reduces  the  risk  of 
overtraining,  while  retaining  nearly  all  of  the  useful  information  in  the  data.  As  such,  linear  models 
are  used  for  our  current  work,  and  rather  than  using  GFA  to  select  what  sensors  to  use  (out  of 
4oC5=658008  possible  combinations),  we  simply  choose  the  most  significant  A  FLD  dimensions, 
where  X  is  allowed  to  vary  between  3  and  1 1 ,  after  which  the  FLD  components  contain  virtually  no 
information  (less  than  1%  of  that  contained  in  the  first  component).  The  value  X  chosen  is  that 


which  minimizes  the  root  mean  squared  error  (RMSE)  of  the  validation  set.  This  process  is  repeated 
for  each  of  the  properties  investigated,  with  a  single  model  chosen  for  each  property,  and  these 
models  are  then  used  to  predict  the  properties  of  those  analytes  in  the  test  set. 

Results 

As  might  be  expected,  some  properties  are  more  easily  related  to  vapor/polymer  interactions,  others 
less.  Results  of  this  analysis  are  shown  in  Table  1.  FCs  refers  to  the  number  of  FLD  components 
used  for  the  model,  and  RMSE/range  refers  to  a  ratio  of  the  RMSE  of  the  test  set  data  to  the  range  of 
the  data.  R2,  slope,  and  int/range  of  the  test  set  all  refer  to  statistics  derived  from  a  correlation  of 
actual  and  predicted  properties  from  the  test  set  data.  Int/Range  refers  to  the  fit  intercept  divided  by 
the  range  of  the  data. 


Property 

#  FCs 

RMSE/Range 

R\  test 

Slope 

Int./Range 

5 

0.0865 

0.864 

1.001 

7.1xl0'4 

van  der  Waals  Volume 

11 

0.1583 

0.607 

0.7750 

0.2489 

van  der  Waals  Area 

10 

0.1607 

0.548 

0.9180 

0.1012 

Radius  of  Gyration 

10 

0.1944 

0.352 

0.6941 

0.3123 

Solubility  Parameter 

9 

0.1116 

0.869 

0.8897 

0.1135 

Hansen  Dipole 

11 

0.1894 

0.513 

0.8037 

0.6141 

Hansen  Electrostatic 

5 

0.2473 

0.626 

1.011 

-0.1589 

Table  1.  Analysis  of  LinearModels  Developed  for  Prediction  of  Thermodynamic/Physical  Parameters 


A  wide  disparity  exists  in  these  results,  with  Kow  and  SP  predicting  very  well,  radius  of  gyration 
predicting  very  poorly,  and  the  other  results  in  between.  Van  der  Waals  area  and  volume  predicted 
well  save  for  an  outlier  which  skewed  the  results.  Note  that  this  method  is  a  fully  blind  experiment 
as  well,  with  the  algorithm  having  never  been  trained  on  the  analytes  whose  properties  it  was  asked 
to  predict,  and  the  test  set  was  twice  as  large  as  the  training  set.  Using  larger  datasets,  or  a  leave- 
one-out  cross  validation  method,9  would  likely  allow  for  greater  success. 
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Section  I 

Response  of  Electronic  Noses  and  the  Analyte  Partition  Coefficient 

S.T.  Lin,  M.  Blanco,  W.  A.  Goddard  III 

Materials  and  Process  Simulation  Center  (MSC),  Caltech 

1.0  Introduction 


It  is  known  that  for  each  polymer  film  sensor  the  response  (increase  in  resistivity 
Ai? )  is  only  a  function  of  the  swelling  volume  of  the  film.1  In  other  words,  the  change  of 
conductivity  is  the  same  as  long  as  the  same  volume  change  AV  in  the  polymer  film  is 
caused,  regardless  of  the  nature  of  the  chemical  molecule. 
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(1) 


This  experimental  evidence  indicates  that  one  could  predict  the  response  of  such 
electronic  nose  if  the  partition  coefficient  of  the  analyte  in  the  polymer  film  and  the  molar 
volume  of  the  analyte  are  known.  Assuming  the  swelling  of  the  polymer  is  ideal,  the 
relative  volume  change  of  the  polymer  film  can  be  determined  from 
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where  Pa  is  the  partial  pressure  of  the  analyte,  Pavap  is  the  vapor  pressure,  V_La  is  the  molar 
volume  of  the  analyte  in  the  pure  liquid  state,  and  the  partition  coefficient  KH  is  defined 
as  the  ratio  of  analyte  concentration  in  the  vapor  phase  c'aap  versus  the  polymer  matrix 
polymer  -n  ^  zero  anaiyte  partial  pressure,  i.e., 

^  polymer 

(3) 
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1.1  Perturbation  Method  for  Calculation  of  Partition  Coefficient 

The  partition  of  analyte  molecules  in  the  gas  and  polymer  is  determined  by  the 
equality  of  the  chemical  potential  in  the  two  phases,  i.e., 

, uvaap(T,P,clap)  =  fj, p°!ymer (T,P, cpolymer ) .  The  Widom’s  particle  insertion  method2  is  one  of 

the  commonly  used  techniques  for  determining  chemical  potentials.  For  a  system 
containing  N&  analyte  molecules  and  Np  polymer  molecules,  the  u  of  analyte  can  be 
expressed  in  terms  of  the  interaction  energy,  ,  between  a  ghost  analyte  and  the 
analyte-polymer  mixture 

fi.(N„Np,P,T )  -  kT In - 2s -  kT\n<—-L - exp(-a'P)  (4) 

<P.  >»..»,  A,  <V 

where  V  is  the  volume  of  the  polymer-analyte  mixture,  q&  and  pa  are  the  internal  partition 
function  and  number  density  of  the  analyte,  Na  and  Np  are  the  number  of  analyte  and 


polymer  molecules,  W  is  the  energy  differences  between  the  systems  with  Na+ 1  analyte 
molecules  and  with  Na  molecules,  i.e.,  T*  =  U(Na  +  l,Np)  -U(Na,Np) ,  and  <>  indicates 

ensemble  averages.  Assuming  the  analyte  is  ideal  in  the  gas  phase,  i.e., 


jiif7  (T,  P)  =  -kT  In — partition  coefficient  can  be  calculated  from  the 


ensemble  average  of  W 
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Molecular  dynamics  simulation  is  a  useful  tool  for  determining  Ku  based  on  the 
Widom’s  formulation  (eq.  5).  Since  preferred  interactions,  i.e.  negative  i|>,  is  rarely 
observed  in  high  density  simulations,  a  common  practice  to  improve  the  ensemble 
average  is  to  gradually  “turn  on”  the  interactions  between  the  ghost  analyte  molecule  and 
the  polymer,  i.e., 


In  Kh 
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where  X  is  a  perturbation  parameter  that  varies  from  0  at  j= 0  to  1  at  j=M,  and  the 
perturbed  interaction  W*  is 

^  =  U(NX  +  -6m,N2)  -  U(NX  +  £j,N2)  =  (e,+1  -  ey.)A Ua_p  (7) 


with  A Ua_p  being  the  interaction  between  the  ghost  (or  perturbed)  analyte  molecule  and 
the  other  molecules  in  the  system,  i.e., 
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where  N\  is  the  total  number  of  atoms  in  the  ghost  analyte  and  N2  is  the  atoms  of  all  other 
atoms  in  the  system.  U^mv  and  UyOUL  are  the  normal  van  der  Waals  and  coulomb 


interactions  between  atoms  i  and  j.  This  is  referred  as  the  Multiple  Step  Thermodynamic 
Perturbation  (MSTP)  method.3  The  MSTP  method  has  also  been  implemented  in 
molecular  dynamics  methods.4 


2.0  Implementation  and  Test  of  the  Multiple  Step  Thermodynamic  Perturbation 
method 

The  MSTP  method  has  been  implemented  in  LAMMPS,  a  parallel  molecular 
dynamics  simulations  code  developed  in  the  Sandia  National  Lab.5  The  following 
analytical  equation  is  used  to  determine  the  values  of  the  perturbation  parameter  X  at  any 
simulation  step  n\ 


e,.  =  e0  exp 
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(9) 


where  no  and  nt  are  the  initial  and  final  simulation  step  number,  and  the  constant  b  is 


determined  by  the  initial  and  final  values  of  X,  i.e,  6  =  1  +  ln-^-.  Equation  9  gives  a 


linear  dependence  of  X  with  n  when  X  is  close  to  1  and  exponential  dependence  when  X  is 
close  to  0.  To  validate  the  LAMMPS  implementation  of  the  MSTP  method,  we  determine 
the  free  energy  of  a  Lennard-Jones  liquid  at  reduced  temperature  of  0.9  and  reduced 
density  of  0.85.  Figure  1  shows  the  instantaneous  and  integrated  values  of  the  free  energy 
of  growing  ( X  varies  from  0  to  1)  and  removing  ( X  varies  from  1  to  0)  a  particle  in  the 
system.  The  results  -2.22  and  2.10  are  in  excellent  agreement  with  the  theoretical  value  of 
±2.14. 


X 


Figure  1.  The  instantaneous  (a)  and  integrated  (b)  values  of  free  energy  change  for  a 
Lennard-Jones  liquid  at  reduced  temperature  of  0.9  and  reduced  density  of  0.85. 


The  free  energy  changes  for  both  growing  (X  varies  from  0  to  1)  and  reducing  (X 
varies  from  1  to  0)  a  ghost  particle  are  shown. 

2.1  Application  to  Electronic  Nose  Sensors:  Partition  Coefficient  of  Analytes  in 
Polymer  Poly-Ethylene  Oxide  (PEO) 

We  aim  at  predicting  the  response  of  electronic  nose  polymeric  sensors  to  a  wide 
variety  of  analytes.  Experimental  data  for  various  polymers,  including  PEO  is  available.1 
To  study  the  partition  coefficient  Ku  of  analyte  in  polymers,  we  performed  MD 
simulations  on  6  analytes  for  which  experimental  data  is  available  (dichloromethane 
CH2CI2,  chloroform  CHCI3,  dibromomethane  CE^B^,  bromoform  CHBr3,  benzene  CgHe, 
and  hexafluorobenzene  in  PEO.  Figure  2  shows  the  typical  change  in  free  energy  as 
a  function  of  X  for  CH2CI2.  (similar  curves  are  observed  for  other  analytes.)  The  free 
energy  curve  for  the  growing  process  starts  at  negative  values,  at  small  values  of  X,  and 
turns  over  to  positive  values,  at  large  X.  This  indicates  a  competition  between  the 
favorable  interactions  and  unfavorable  repulsive  interactions  at  different  values  of  X.  At 
small  X  values  the  interaction  between  the  analyte  and  polymer  is  weak.  The  fractional 
analyte  can  occupy  the  free  volume  of  the  polymer  without  changing  the  polymer 
conformation.  The  increase  of  X  improves  the  attractive  interactions  between  the  analyte 
and  polymer,  leading  to  negative  values  of  the  free  energy.  However,  at  large  values  of  X, 
the  increasing  repulsive  interactions  become  more  significant.  In  other  words,  the 
polymer  has  to  change  its  conformation  to  create  more  volume  for  the  growing  analyte. 
Therefore,  a  positive  free  energy  is  observed  at  large  X  values. 


Figure  2.  Instantaneous  values  of  free  energy  for  dichloromethane  in  PEO  at 
different  X. 


The  predicted  partition  coefficients  for  all  six  analytes  are  compared  with 
experiment  in  Table  1.  A  qualitative  agreement  between  experimental  measurements  is 
obtained.  Discrepancies  may  be  the  result  of  insufficient  sampling,  which  is  known  to  be 
crucial  for  polymer  simulations.  At  the  same  time  the  measurements  are  subject  to 
experimental  uncertainty.  For  example,  ellipsometry  measurements  are  typically  higher, 
by  factors  of  3,  than  the  results  of  the  quartz-crystal  microbalance  (QCM)  quoted  here. 


Table  1.  Comparison  of  the  predicted  and  experimental  values  of  partition 


coefficien 


ln/fn  for  6  analytes  in  PEO  at  300  K 


analyte 

EE3S3SS3SHMH 

CH2CI2 

1.91 

4.85 

CHCI3 

3.34 

6.44 

CH2Br2 

4.47 

6.36 

CHBr3 

6.42 

7.39 

c6h6 

2.58 

5.89 

c6f6 

0.80 

-2.01 

To  examine  the  importance  of  sampling,  we  have  also  performed  two  tests:  (1) 
increase  the  number  of  simulation  cycles  and  (2)  increase  the  number  of  perturbed 
analyte  molecules  in  the  system.  Table  3  compares  the  results  using  either  one  or  three 
cycles,  i.e.,  repeating  the  change  of  X  from  0  to  1  and  back  to  0  one  or  three  times.  It  is 
found  that  the  increase  of  simulation  cycles  uniformly  decreases  the  values  and  improves 
the  prediction  of  lnifn-  Therefore,  the  increase  of  cycles  can  reduce  the  bias  towards  poor 
polymer-analyte  interactions.  It  seems  then  that  the  polymer  undergoes  a  non-reversible 
hysteresis  that  requires  a  few  cycles  to  equilibrate. 

Table  2.  The  predicted  partition  coefficient  IhA'h  with  multiple  MSTP  cycles. 


analyte 

lnXnCexpt) 

lnKntcalc,  1  cycle) 

lnifn(calc,  3  cycles) 

CH2CI2 

1.91 

4.85 

0.96 

CH2Br2 

4.47 

6.36 

3.83 

CHBr3 

6.42 

7.39 

4.75 

A  second  test  is  performed  on  the  CH2CI2-PEO  system  using  1,5,  and  10  analyte 
molecules.  The  multiple  analytes  behave  as  individual  probes  sampling  different  regimes 
of  the  polymer.  Table  3  compares  the  predictions  of  IilKh  using  different  number  of 
analyte  probes.  It  is  found  that  having  multiple  analytes  in  the  system  also  help  improve 
the  sampling.  This  approach  is  much  more  efficient  than  the  previous  multiple  cycle 
approach  and  is  preferred  for  the  future  study  of  the  analyte  partition  in  polymer  sensors. 


Table  3.  Comparison  Effects  of  samplings  and  number  of  analy 


analyte 

CH2CI2 

lnKH(expt) 

1.91 

1ilKh(1  CH2C12,  1  cycle) 

4.85 

ln£H(l  CH2CI2,  3  cycles) 

0.96 

ln£H(5  CH2CI2,  1  cycle) 

0.80 

lni^lO  CH2CI2,  1  cycle) 

1.12 

es  in  the  simulation. 


Conclusions 

In  this  work  we  have  used  the  thermodynamic  perturbation  method  to  determine 
the  partition  coefficient,  the  key  component  of  determining  the  polymer  sensor  response, 
for  six  analytes  in  poly  (ethylene  oxide).  We  have  implemented  the  multiple  step 
thermodynamic  perturbation  (MSTP)  algorithm  in  a  parallel  molecular  dynamic 
simulation  program  LAMMPS.  In  MSTP,  we  adjust  the  polymer-analyte  interaction 
through  a  perturbation  parameter  X  whose  value  varies  between  0  and  1 .  Integrating  the 
free  energy  change  of  the  polymer-analyte  composition  from  X=0  to  1  gives  the  partition 
coefficient.  Good  agreement  with  experimental  measurements  is  found  when  enough 
samplings  are  considered.  The  MSTP  method  is  a  practical  and  useful  tool  for  studying 
the  dependence  of  polymer  sensor  repose  on  different  analyte  and/or  polymer  chemical 
structures,  humidity,  temperature,  and  plasticitizers. 
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Section  II 

Hildebrand  and  Hansen  Solubility  Parameters  from  Molecular  Dynamics  with 
Applications  to  Electronic  Nose  Polymer  Sensors 

M.  Blanco,  and  W.  A.  Goddard  III 

Materials  and  Process  Simulation  Center,  California  Institute  of  Technology 
1.0  Introduction 

In  1936  Joel  H.  Hildebrand  proposed1  a  simple  definition  for  a  “solubility  parameter”  that  would 
provide  a  systemic  description  of  the  miscibility  behavior  of  solvents  and  which  subsequently  has 
found  multiple  uses  in  chemistry.  This  solubility  parameter  6  is  defined  as  square  root  of  the 
cohesive  energy  density,  the  heat  of  vaporization  divided  by  the  molar  volume.  Hansen2  proposed 
an  extension  of  the  Hildebrand  parameter  method  to  estimate  the  relative  miscibility  of  polar  and 
hydrogen  bonding  systems.  In  Hansen’s  approach  the  Hildebrand  solubility  parameter  is  split 
into  three  components,  polar,  dispersion,  and  hydrogen  bonding,  thus  the  name  3D  solubility 
parameters.  The  three  components  are  empirically  fitted  to  define  the  miscibility  characteristics 
of  the  solvent.  Solvents  with  similar  Hansen  solubilities  are  miscible  in  most  proportions; 
dissimilar  values  yield  limited  solubilities.  Hildebrand  and  Hansen  solubility  parameters  are 
useful  for  selecting  solvents  and  additives  in  formulations,  for  the  blending  of  polymers,  for  the 
control  of  kinetics  and  monomer  sequence  distributions  in  copolymers,  and  for  the  proper 
selection  of  time-release  formulations  in  the  delivery  of  pharmaceuticals.  The  first  principles 
predictions  of  Hildebrand  solubility  parameters  is  of  great  practical  value  in  chemical  formulation 
work.  Unfortunately,  it  has  been  difficult  to  obtain  sufficient  experimental  information  to  define 
all  of  the  parameters  needed  in  developing  new  formulations. 

Consequently,  several  groups  have  attempted  to  develop  general  computational  tools  to  estimate 
Hildebrand  and  Hansen  solubility  parameters.  Choi  and  Kavasallis  first  used  atomistic 
simulations  to  estimate  the  solubility  parameters  of  a  class  of  alkyl  phenol  ethoxylates3  and  later 
applied  it  to  the  estimation  of  3D  Hansen  solubility  parameters4.  A  related  method  has  been 
applied  to  the  estimation  of  the  solubility  parameters  for  distributions  of  asphaltenes,  resins,  and 
oils  from  crude  oils  and  related  materials.5  The  accuracy  of  these  methods  depends  on  the  correct 
building  of  the  bulk  structure  as  well  as  on  the  molecular  force  field  parameters  used  in  the 
calculations.  Numerous  approaches  for  building  amorphous  polymers  and  liquids  have  been 
published.6'11  Some  of  these  methods  involve  growing  the  polymer  chains  at  a  fixed 
experimental  density  using  rotational  isomeric  state  (RIS)  statistics  in  combination  with  a  scaled 
down  atomic  radius  followed  by  potential  energy  minimization  with  periodic  boundary 
conditions.  Other  methods  simulate  a  “polymerization”  process  to  grow  the  chain  in  a  fixed 
density.  A  computationally  expensive  protocol  involving  chain  growth  at  low  density  followed 
by  a  pressure- induced  compression  with  molecular  dynamics  has  also  been  reported.12  Most  of 
these  methods  have  been  successfully  used  to  generate  amorphous  structures  and  have  correctly 
predicted  the  solubility  parameters  of  a  few  polymers. 

Here  we  report  on  a  multi-sample  molecular  dynamics  method  which  provides  a  feasible  tool  for 


estimating  Hildebrand  and  Hansen  solubility  parameters  without  the  need  for  experimental  data. 
The  molecular  dynamics  method  developed  in  this  work  is  particularly  useful  in  rapidly 
generating  structures  of  polymers  with  large  monomer  units  containing  rings  or  other  complex 
groups.  The  finite  number  of  densification  and  equilibration  steps  regardless  of  polymer  size 
allows  for  a  gradual  packing  adjustment  and  the  uniform  redistribution  of  stresses  among  the 
polymer  chains.  This  new  method  was  validated  by  several  studies  where  solubility  parameter 
calculations  were  successfully  correlated  with  experimental  measurements. 

For  improved  accuracy,  the  new  method  employs  quantum  mechanical  charges  of  single 
molecules.  However,  semi-empirical  methods  for  charge  assignment,  such  as  QEq,  give 
somewhat  comparable  results  for  molecules  containing  first  group  elements.  The  most 
significant  approximation  comes  from  the  use  of  a  generic  force  field  for  the  estimation  of 
dispersion  and  hydrogen  bonding.  Approximations  not  withstanding,  calculated  Hildebrand 
parameters  compare  well  with  experimental  values  for  a  series  of  solvents  and  monomer 
molecules.  As  an  application  example,  we  illustrate  the  use  of  these  values  in  the  design  of 
polymer  sensors  for  an  electronic  nose. 


1.1  Method 

The  Hildebrand  solubility  parameter  for  a  pure  liquid  substance  is  defined  as  the  square  root  of 
the  cohesive  energy  density. 

6  =[(AHv-RT)/Vm]m  (1) 

AHV  is  the  heat  of  vaporization  and  Vm  the  molar  volume.  Typical  units  are 

1  hildebrand  =  1  cal 12 crn3/2  =  0.48888  x  MPa  in=  2.4542  xlO  2  (Kcal/mol) 1/2  A'm  (la) 


Hansen  proposed  an  extension  of  the  Hildebrand  parameter  to  estimate  the  relative  miscibility  of 
polar  and  hydrogen  bonding  systems 

<52  =  V  +  V  +  <5* 2  (2) 

where  $  ,  $ and  6h  are  the  dispersion,  electrostatic,  and  hydrogen  bond  components  of 
6,  respectively.  For  molecules  whose  heats  of  vaporization  can  be  measured  or  calculated,  one 
can  easily  determine  the  value  of  5.  The  Hansen  solubility  parameters  in  Equation  (2)  are 
determined  empirically  based  on  multiple  experimental  solubility  observations.  However,  for 
polymers  the  Hansen  parameters  are  assigned  to  the  parameters  of  the  solvent  causing  the 
maximum  swelling  in  a  series  of  polymer  swelling  experiments.  Thus,  the  two  quantities 
represented  by  Equation  (1)  and  (2)  are  expected  to  be  similar  but  not  identical,  because  the 
Hildebrand  parameters  are  not  always  determined  from  heats  of  vaporization,  particularly  for 
substances  with  high  boiling  points.  For  polymers,  a  variety  of  other  experimental  methods  are 
also  used14  leading  to  a  wide  range  of  values  for  any  given  polymer. 

We  report  here  the  Cohesive  Energy  Density  (CED)  protocol  based  on  Molecular  Dynamics 
(MD)  calculations  under  periodic  boundary  conditions  (PBC)  to  determine  an  ensemble  of 


temperature  and  pressure  equilibrated  structures  from  which  we  can  extract  properties  of  the 
liquid,  inlcuding  the  Hansen  and  Hildebrand  solubilities  of  solvents  and  polymers.  CED  leads  to 
sytematic  estimates  of  the  uncertainties  in  these  quantities  (within  the  model  and  the  size  of  the 
ensemple)  that  are  often  better  than  the  experimental  ones.  .  The  CED  method  overcomes  the 
common  equilibration  problems  with  condensed  phase  molecular  dynamics,  i.e.,  how  to  choose 
initial  molecular  configurations  not  far  from  equilibrium  at  normal  densities.  Significant  amounts 
of  simulation  time  are  usually  required  to  equilibrate  the  badly  packed  molecules  often  generated 
with  Monte  Carlo  methods.  In  particular  to  obtain  densely  packed  polymers  with  their 
enormosous  number  of  torsional  degrees  of  freedom  ,  often  leads  to  highly  nonequilbrium 
dihedral  populations.  Thus,  care  must  be  taken  to  generate  an  ensemble  of  thermally  accessible 
conformations  not  far  from  equilibrium.  These  two  requirements,  condensed  phase  densities  and 
equilibrated  molecular  conformations,  are  satisfied  through  the  following  method15 

1.  A  cubic  periodic  unit  cell  containing  a  given  number  of  molecules  is  built  at  50%  of  the 
target  density.  Generally  four  polymer  chains  are  sufficient  although  for  very  high 
molecular  weights  even  one  chain  can  be  adequate.  For  the  solvents,  between  16  to  64 
solvent  molecules  are  adequate.  We  find  that  for  building  the  structure  it  is  useful  to  scale 
Van  der  Waals  radii  by  a  factor  of  0.30  to  get  initial  strucures  that  will  eventually  lead  to  a 
good  ensemble.16  For  convenience  in  this  paper,  we  used  the  experimental  density  of  the 
solvents  and  polymers  as  a  target  value  since  they  are  available,  however,  systems  with 
unknown  density,  would  only  require  one  additional  step  of  NPT  dynamics  to  obtain  the 
trial  density  for  CED. 

2.  The  initial  polymer  amorphous  structures  are  constructed  using  the  rotational  isomeric 
state  (RIS)  table  and  a  suitable  Monte  Carlo  procedure16  to  achieve  a  correct  distribution 
of  conformational  states  in  the  low  density  sample. 

3.  The  charges  of  the  isolated  solvent  or  polymer  molecules  are  defined  using  the  charge 
equilibration  method13  or  are  obtained  from  quantum  mechanical  calculations. 

4.  The  force  field  parameters  are  taken  from  a  suitable  force  field,  such  as  the  generic 
Dreiding  forcefield6,  Universal  force  field  (UFF35)  etc 

5.  Minimization:  The  potential  energy  of  the  bulk  system  is  minimized  for  5000  steps  or 
until  the  atom  rms  force  converges  to  0.10  kcal/mol-A. 

6.  Dynamics:  750  steps  of  Molecular  Dynamics  (1  femtosecond/step)  at  a  temperature  of 
700  K  using  canonical  fixed  volume  dynamics  (NVT)  are  carried  out  to  anneal  the 
sample. 

7.  Compression:  The  reduced  cell  coordinates  are  shrunk  such  that  the  density  is  64%  of  the 
target  density. 

8.  The  atomic  coordinates  are  minimized  and  dynamics  is  done  on  the  system  with  the 
previously  described  procedure  holding  the  cell  fixed  (steps  5-6). 

9.  A  total  of  five  compression,  minimization,  and  dynamics  cycles  are  performed  until  the 
density  reaches  120  %  of  the  target  density. 

10.  The  cell  parameters  are  increased  in  five  cycles  of  expansion,  minimization  and  dynamics, 
until  the  target  density  is  reached. 

1 1 .  Finally,  the  sample  is  allowed  to  relax  in  a  minimization  involving  the  cell  and  the  atomic 


coordinates. 


12.  Molecular  dynamics  are  performed  for  20  picoseconds.  The  first  10  picoseconds  are  used 
for  thermalization  of  the  sample  at  the  desired  temperature.  The  last  10  picoseconds  are 
used  for  averaging  of  cell  volume  and  potential  energies  (dispersion  of  van  der  Waals, 
polar  or  Coulomb,  and  hydrogen  bonding). 

13.  The  Hansen  enthalpy  components  are  calculated  by  subtracting  the  potential  energy  of  the 
bulk  system  from  the  sum  of  the  potential  energies  of  the  individual  molecules  as  if 
separated  by  an  infinite  distance. 

14.  This  process  is  repeated  ten  times  with  different  initial  random  conformations  and 
packing. 

15.  Hansen  solubility  parameters  and  molar  volumes  are  computed  as  well  as  their  standard 
deviations.  95%  confidence  limit,  two  standard  deviations  from  the  average  value,  are 
used  to  identify  outliers.  Typically  none  or  at  most  two  outliers  are  encountered  in  a  ten 
sample  run. 

Hildebrand  and  Hansen  solubilities  are  calculated  from  the  molecular  dynamics  average  energy 
components  of  the  condensed  phase  simulation,  single  unit  cell  Ec,  the  energy  components  of  the 
individual  molecules,  and  the  volume  of  the  simulated  sample,  Vc  as  follows 
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where  <  >  indicates  a  time  average  over  the  duration  of  the  dynamics,  n  the  number  of  molecules 
and  N0  is  Avogadro’s  number.  When  E,  is  taken  to  be  the  total  energy,  the  value  of  6  is  the 
Hildebrand  solubility  parameter,  or  one  of  the  components,  E^.^ond,  Ecoulombi  E diSpersiom  which 
results  in  the  Hansen  solubility  parameters.  We  stress  the  importance  of  using  the  thermally 
equilibrated  ensemble  of  molecular  conformations  to  estimate  the  gas  phase  term  £,,  instead  of  n 
times  the  minimized  energy  of  one  molecule.  For  simplicity  the  gas  phase  ensemble  average  is 
approximated  by  the  average  of  the  isolated  molecules  taken  from  the  condensed  phase 
simulation,  averaged  over  the  entire  molecular  dynamics  at  the  desired  temperature. 

The  builder  converts  an  existing  model  into  an  amorphous  structure  by  manipulating  the  model’s 
rotatable  bonds.  Each  unique  torsion  can  be  defined  using  a  Monte  Carlo  procedure  with 
statistical  weights  given  by  a  previously  built  rotational  isomeric  state  table.16  Conformations 
are  rejected  if  two  or  more  atoms  come  in  contact  closer  than  a  van  der  Waals  scale  distance.  The 
resulting  amorphous  structure  can  then  be  relaxed.  In  polymer  calculations,  the  number  of 
monomers  in  each  chain  is  usually  determined  such  that  the  total  volume  of  the  four  chains  is 
approximately  5900  A3.  Alternatively,  a  degree  of  polymerization  of  30  suffices  to  give 
converged  values.  In  such  polymer  samples,  the  minimum  number  of  atoms  is  around  1000.  The 
overall  procedure  is  illustrated  in  Figure  1 . 


2.0  Results 


Over  sixty  common  solvents,  reactants,  and  monomers  of  significantly  different  structure, 
polarity,  and  chemical  composition  were  chosen  from  various  experimental  compilations  of 
Hildebrand  solubility  values  available  in  the  literature.  Quantum  Mechanical  (Hartree  Fock  6- 
31G**  full  geometry  optimization,  Mulliken  and  electrostatic  potential  (ESP)  with  a  constraint  to 
reproduce  the  quantum  dipole  moment  from  the  wavefunction)  charges  were  assigned  to  the 
molecules  for  ten  separate  independent  molecular  dynamics  simulations.  Table  I  (attached  after 
references)  contains  the  averaged  results  of  1 0  CED  simulations  for  the  two  charge  assignment 
methods  for  each  solvent  molecule.  For  comparison  the  various  experimental  values  found  in 
the  literature  and  their  deviations  are  included.  Table  II  gives  an  example  of  the  output  of  the 
CED  procedure. 

Figures  2  and  3  show  the  correlations  between  experimental  and  predicted  values.  Based  on  this 
data  the  accuracy  of  the  CED  method  is  0.8  hildebrands.  The  experimental  RMS  value  calculated 
from  the  various  literature  sources  is  0.43  hildebrands.  Although  most  solvents  fall  within  the 
experimental  error  a  few  predictions  are  clearly  outside  the  range  of  measured  values.  Exclusion 
of  the  six  worst  cases  (formic-acid,  acetic  acid,  dichlorodifluoromethane,  acrylic  acid,  methyl 
formamide,  and  malononitrile)  from  the  predictions  reduces  the  RMS  to  0.6  hildebrands. 

The  accuracy  of  the  molecular  dynamics  results  directly  depends  on  the  accuracy  of  the  intra  and 
inter-molecular  potential  atomic  parameters  (force  field)  and  to  some  extent  on  the  modeling 
protocol.  This  problem  is  in  part  overcome  with  force  fields  that  accurately  reproduce  the 
experimentally  measured  bond  distance,  angles,  as  well  as  the  respective  force  constants  of  small 
molecules.  Less  effort  has  gone  into  optimizing  the  van  der  Waals  parameters  in  such  force 
fields.  Precision  on  the  other  hand  is  strongly  dependent  on  the  molecular  dynamics  procedure 
employed  to  prepare  the  samples.  No  significant  differences  in  precision  were  found  for  the 
worst  six  cases  (0.60,  0.62,  0.40,  0.37,  0.57,  0.42  hildebrands  respectively)  when  compared  to  the 
average  precision  across  all  solvents  (0.44  hildebrands). 

We  speculate  that  the  assigned  van  der  Waals  force  field  parameters  (our  generic  force  field  was 
not  particularly  fitted  to  halogens  and  nitrogen  containing  compounds)  play  a  role  in  the  accuracy 
of  our  predictions.  We  made  no  attempt  to  adjust  the  force  field  parameters  here.  However  we 
point  to  the  possibility  of  using  the  CED  method  together  with  experimental  heats  of  vaporization 
and  densities  for  the  estimation  of  van  der  Waals  Lennard-Jones  parameters  and/or  the  hydrogen 
bond  terms  for  the  various  chemical  atom  types  represented  by  these  compounds.  For  example, 
systematic  underestimation  of  the  solubility  parameter  is  observed  through  the  calculated  versus 
experimental  ratio  6esp/Sexp  for  alcohols  and  amides  (  2-ethyl- lhexanol  0.89,  2-ethyl- 1 -butanol 
0.92,  1-pentanol  0.95,  n-butanol  0.92,  n-propanol  0.91,  furfuryl  alcohol  0.96,  ethanol  0.87,  1,3- 
butanediol  0.91,  methanol  0.89,  N,N-dimethylacrylamide  0.91,  dimethylacetamide  0.94, 
dimethylformamide  0.87,  methylformamide  0.84).  This  suggests  that  the  Dreiding  parameters  for 
the  H-bond  term  (Do=2.5  Kcal/mol,  Ro=3.2  A)  could  be  modified  to  increase  the  accuracy  of  the 
predictions.  Both  parameters  may  be  involved  since  the  density  of  these  two  groups  of 
compounds  is  also  underestimated.  Finally,  there  seems  to  be  a  systematic  overestimation  of  the 
solubility  parameter  for  the  organic  acids.  The  Sesp/^exp  ratio  is  consistently  high  (propionic  acid 
1.18,  acetic  acid  1.30,  methacrylic  acid  1.04,  formic  acid  1.30,  acrylic  acid  1.21).  The  effect  is 


opposite  the  hydrogen  bond  effect  previously  mentioned.  We  assume  that  the  molecules  in  the 
gas  phase  are  non-interacting.  Many  low  molecular  weight  acids  exist  in  a  dimerized  form  in  the 
gas  phase.  If  we  assume  that  half  of  the  intermolecular  hydrogen  bonds  are  preserved  in  the  gas 
phase,  the  solubility  parameter  will  be  decreased  by  about  3  hildebrands,  bringing  theory  and 
experiment  quite  close  in  agreement.  We  now  discuss  the  effect  of  charge  assignment  methods. 

Although  Mulliken  charges  are  quite  useful  in  determining  the  structure  of  molecules  in  the  gas 
phase  these  appear  to  be  less  accurate  than  Electrostatic  Potential  charges  (ESP)  for  the 
determination  of  condensed  phase  properties.  We  advocate  the  use  of  ESP  charges  for  the 
estimation  of  solubility  parameters  and  the  use  of  Mulliken  charges  for  conformational  studies  in 
the  gas  phase.  It  appears  that  the  far  field  representation  of  the  ESP  charges  captures  more 
accurately  the  physical  interactions  between  molecules  in  the  condensed  phase. 

We  compare  our  Molecular  Dynamics  results  to  other  predictive  methods  available  in  commercial 
software  packages  such  as  Synthia-Fedors  and  Synthia-van-Krevelen.25  These  methods  can  be 
considered  state-of-the-art  group  additivity  methods,  relying  on  topological  descriptors  and  other 
single  molecule  quantities  to  make  predictions  based  on  correlations  and  parameter  extractions 
from  large  databases  of  solubility  parameters.  Although  intended  for  predictions  on  polymers 
these  parametric  methods  require  only  the  isolated  monomer  structure  to  make  a  prediction.  The 
methods  are  fast  and  simple  to  use.  In  contrast,  the  MD  method  presented  here  requires  a  full 
condensed  phase  simulation  of  the  compound  of  interest.  Nonetheless  the  CED  Molecular 
Dynamics  method  is  non-parametric.  Beyond  the  predetermined  force  field,  in  our  case  a  generic 
force  field  published  in  the  mid  80’s,  CED  used  no  adjustable  parameters  and  no  experimental 
input  information.  Moreover,  in  principle  the  CED  molecular  dynamics  method  can  make 
predictions  as  a  function  of  pressure  and  temperature  and  it  is  general  enough  to  deal  with 
complex  mixtures,  including  solvent/polymer  mixtures.  The  average  root  mean  square  (RMS) 
difference  for  the  CED  results  and  for  two  group  additivity  predictive  methods,  Synthia-Fedors 
and  Synthia-van-Krevelen5  with  respect  to  the  experimental  values  is  shown  Table  III 

Solubility  parameters  predicted  employing  the  current  molecular  dynamics  simulation 
methodology  can  also  be  compared  to  those  calculated  employing  molecular  dynamics 
simulations  by  Rigby  et  al.26a.  Employing  the  PCFF  forcefield  and  Amorphous  Cell/Discover 
programs,25  Rigby  et  al.  predict  solubility  parameters  for  13  of  the  molecules  in  Table  I  which 
have  an  average  absolute  difference  with  experiment  of  0.92  (cal/cm3)1 2  which  is  similar  to  that 
predicted  by  the  current  methodology  (absolute  difference  of  1.1  for  the  13  molecule  set).  As 
with  the  current  simulations,  the  largest  differences  in  the  set  are  seen  for  two  acid  molecules. 
The  current  authors  have  calculated  the  solubility  parameters  for  hexane,  acetone,  and  n-propanol 
employing  molecular  dynamics  simulations  with  the  COMPASS  forcefield  and  Amorphous 
Cell/Discover  programs25'  Average  absolute  differences  with  experiment  for  this  set  of  three  was 
found  to  be  0.25  with  COMPASS  methodology  which  is  slightly  smaller  than  the  0.55  difference 
observed  for  the  current  CED  methodology.  The  COMPASS  forcefield  has  been  extensively 
optimized  to  reproduce  the  heats  of  vaporization  of  a  large  number  of  organic  liquids.  For 
example,  in  a  related  COMPASS  forcefield  study, 26c  Sun  has  calculated  heats  of  vaporization  for 
100  compounds  to  within  an  average  percent  error  of  experiment  of  -0.2%  with  maximum  errors 
of  14.6%  and  14.5%.  This  is  to  be  contrast  with  our  current  approach  where  generic  force  field 
was  employed  without  any  further  optimization.  Eichinger  and  coworkers26b  have  employed  the 
COMPASS  forcefield  to  compute  solubility  parameters  for  Ultem  oligomers,  related  molecules, 
and  solvent  molecules  including  toluene.  Not  surprisingly  the  calculated  solubility  parameter  for 
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toluene  is  0.07  (cal/cm  )  closer  to  the  average  experimental  value  (8.94  (cal/cm  )  )  than  the 
previous  PCFF  forcefield  value.263  The  calculated  toluene  solubility  parameter  value  of  9.0 
(cal/cm  )  compares  well  to  the  current  ESP  calculated  solubilty  parameter  of  9.23  (cal/cm  )  . 


2.1  Application  to  the  Electronic  Nose 


An  electronic  nose  has  been  built  at  Caltech27'29  employing  an  array  of  polymer  sensors.  Sensors 
are  built  with  conducting  leads  connected  through  thin  film  polymers  loaded  with  carbon  black. 
Odorant  detection  relies  on  a  change  in  electric  resistivity,  AR/R,  of  the  polymer  film  as  function 
of  the  amount  of  swelling  caused  by  the  odorant  compound.  The  amount  of  swelling  depends 
upon  the  chemical  composition  of  the  polymer  and  the  odorant  molecule.  An  array  of  twenty 
carbon  black  loaded  polymers  give  rise  to  a  specific  change  in  resistivity  patterns  upon  exposure 
to  a  given  molecular  species.  The  pattern  is  unique  and  unambiguously  identifies  the 
compound  28,29 

The  experimentally  determined  changes  in  relative  resistivity,  AR/R,  of  seven  polymer  sensors 
upon  exposure  to  twenty-four  solvent  vapors  was  correlated  with  the  calculated  Hansen  solubility 
components.  The  permeability  of  a  given  odorant  in  a  polymer  is  given  by20 


P  =  Aexp 


A  H, 

kBT 


(4) 


where  A  is  the  pre-exponential  factor  related  to  entropy,  AHS  is  the  heat  of  sorption  of  the  solute 
and  Ed  is  the  activation  energy  for  diffusion  of  the  molecule  in  the  polymer.  We  assume  that  the 
relative  change  in  resistivity  is  directly  proportional  to  the  odorant’s  permeability.  The  following 
expression  was  used  to  correlate  AR/R  with  the  Hansen  components  of  the  cohesive  energy  of  the 
polymer  and  solvent  as  well  as  the  molar  volume  of  the  solvent. 

AR/S.R,  exp(-rr,)exp[V  AW  -  «,')]  (  5  ) 


yVs  is  the  activation  energy  of  diffusion  of  the  solute  in  the  polymer,  proportional  to  the  molar 
volume  of  the  odorant,  Vs.  The  exponential  factor  y  is  a  best-fit  parameter.  We  base  this  relation 
on  the  experimental  observation  that  the  diffusion  coefficients  of  various  molecules  is  linearly 
related  to  the  molar  volume  of  the  solute  in  the  case  where  the  actual  temperature  is  greater  than 
the  glass  transition  (Tg)  of  the  polymer.31  This  approximation  is  used  in  our  analysis  regardless  of 
Tg.  df  (i=l,2,3)  are  the  cohesive  energy  density  component  of  the  solvent  s,  where  i=l,  2  and  3 
refer  to  the  electrostatic,  dispersion  and  hydrogen  bond  components  respectively.  Similarly,  <5,p  is 
the  i-th  cohesive  energy  component  of  the  polymer  sensor  p.  The  exponential  coefficients  /3,  are 
treated  as  best  fit  parameters  as  well  as  is  the  pre-exponential  term  Ro-  It  should  be  noted  that  we 
preserve  the  sign  of  the  energy  components  in  equation  (5),  usually  lost  in  the  definition  of 
Hansen  and  Hildebrand  parameters.  This  is  important  because  such  interactions  can  be  attractive 
or  repulsive,  depending  on  the  polymer/odorant  mixture  in  question.  For  simplicity  we  employed 
the  charge  equilibration,  Qeq,13  method  to  assign  atomic  charges  to  polymers  and  solvent 
molecules. 


The  results  of  fitting  equation  (5)  to  experimental  changes  in  resistivity  are  shown  in  Figure  4 
for  seven  electronic  nose  polymer  sensors  and  twenty-four  solvents.27  Pearson’s  correlations 
between  the  experimentally  determined  change  in  resistivity  and  the  Hansen  solubilities  are 
shown  for  polymer  sensors  (poly(methylmethacrylate)  (PMMA),  poly(4-hydroxystyrene) 
(P4HS),  polyethyleneoxide  (PEO),  polyethylene  (PE),  poly(ethylenevinyl  acetate)  (PEVA), 
polysulfone,  and  caprolactone)  in  Table  IV.  The  calculated  Hansen  solubilities  for  the  seven 
polymers  and  twenty-four  solvents  are  summarized  in  Table  V. 

The  correlation  was  particularly  good  for  polysulfone,  poly(4-hydroxystyrene)  and  PEVA 
(polyethylene-co-vinyl  acetate)  and  especially  poor  for  polymethylmethacrylate  based  on  both 
correlation  slope  and  the  Pearson  R  values  for  the  linear  fit.  Polysulfone  appears  to  discriminate 
between  solvents  of  different  sizes  since  the  free  volume  fraction  is  small  and  the  free  volume 
distribution  may  be  narrow,  resulting  in  a  "molecular  sieve"  effect.  Additionally,  the 
experimental  relative  change  in  resistivity  in  polysulfone  ranges  from  zero  to  1.0,  which  makes  it 
a  particularly  good  high-resolution  sensor. 

The  polyethylene-co-vinyl  acetate  detector  also  correlates  reasonably  well  with  the  theoretical 
relative  change  in  resistivity.  However,  the  relative  change  in  resistivity  range  is  smaller 
compared  to  polysulfone  indicating  that  it  is  less  discriminating  towards  ester  and  alcohol 
solvents.  A  possible  explanation  that  accounts  for  this  observation  is  that  PEVA  contains  polar 
ester  functional  groups  due  to  the  vinyl  content  (18  %),  as  well  as  non-polar  components  due  to 
the  polyethylene  content  (82  %).  PEVA  has  a  glass  transition  below  room  temperature  and  as  a 
result  contains  a  large  free  volume  fraction.  This  decreases  the  sensitivity  towards  molecules  of 
different  sizes  compared  to  high  Tg  polymers  such  as  polysulfone.  The  third  particularly  good 
detector  in  terms  of  signal  correlation  with  theoretical  prediction  is  poly(4-hydroxy  styrene). 

This  detector  is  particularly  sensitive  to  molecules  functionalized  with  highly  polar  groups  such 
as  alcohol  due  to  the  hydroxyl  functional  group.  However,  the  sensitivity  of  this  sensor  to 
moderately  polar  or  non-polar  solvents  such  as  esters  is  particularly  low. 

3.0  Conclusions 

The  CED  Molecular  Dynamics  method  presented  here  provides  a  reproducible  method  to  predict 
the  Hildebrand  solubility  parameters  for  polymers  and  solvents.  Its  precision  (0.44  hildebrands) 
is  comparable  to  the  experimental  error  (0.43).  We  investigated  the  use  of  compression  and 
expansion  cycles,  simulated  annealing,  charge  assignment  methods,  and  statistical  sample 
averaging.  It  is  important  to  start  from  a  low-density  sample  to  achieve  equilibrated 
conformational  statistics  within  a  reasonable  computational  time.  Ten  samples,  with  roughly 
1,000  atoms  each,  seem  adequate  to  estimate  these  properties. 

The  results  compare  well  with  experimental  values,  although  accuracy  (0.88  hildebrands)  is 
somewhat  lower  than  precision  (0.44  hildebrands),  mostly  a  deficiency  of  the  force  field.  No 
attempt  was  made  to  refine  the  force  field  parameters  to  improve  the  accuracy  of  the  method, 
although  such  a  possibility  is  clearly  present.  The  present  method  should  be  added  to  the 
existing  methods  (ab  initio  EOS,  x-ray  structure  fittings)  currently  used  to  estimate  Lennard- 
Jones  parameters  or  hydrogen  bond  parameters. 


The  CED  Molecular  Dynamics  calculated  Hildebrand  and  Hansen  parameters  give  systematic 
measures  of  the  solvency  of  a  given  chemical  compound.  Such  values  are  of  practical  use  for  the 
design  of  new  materials,  formulations,  and  processes.  As  parallel  molecular  dynamics 
algorithms  are  implemented  simulation  times  will  be  reduced  and  the  prediction  of  solubility 
parameters  will  become  even  more  practical.  For  example,  the  estimation  of  Hildebrand  and 
Hansen  solubility  parameters  takes  approximately  two  hours  in  a  dual  processor  Linux  computer. 
Using  a  highly  parallel  particle-mesh  algorithm32'33  to  integrate  the  dynamics,  the  CPU  times 
can  be  reduced  to  ten  minutes  using  twenty  four  processors.  In  such  a  computational 
environment,  it  becomes  practical  to  automate  the  population  of  databases  of  solvents  and 
complex  mixtures  with  Hildebrand  and  Hansen  solubility  parameters  for  product  formulation 
work. 


3.1  Acknowledgements 

This  work  was  supported  by  an  ARO/MURI  and  in  part  by  a  grant  by  the  3M  Company  and 
Owens  Corning.  The  facilities  of  the  MSC  were  partly  funded  by  NSF  MRI  and  ARO/DURIP 
and  are  also  supported  by  grants  from  DOE- ASCI,  Chevron,  NIH,  ONR,  Seiko-Epson,  Avery- 
Dennison,  Kellogg’s,  General  Motors,  Beckman  Institute,  Asahi  Chemical,  and  Nippon  Steel. 


References 


1.  J.H.  Hildebrand,  The  Solubility  of  Non-Electrolytes.  New  York:  Reinhold,  (1936). 

2.  Hansen,  C.M.,  J  Paint  Technol  39,  511,  (1967). 

3.  P.  Choi,  T.A.  Kavassalis,  and  A.  Rudin,  J.  Coll.  Interface  Sci.,  150,  386  (1992). 

4.  T.A.  Kavassalis,  P.  Choi,  and  A.  Rudin,  Molecular  Simulation  11,  229  (1993). 

5.  E.  Rogel,  Energy  &  Fuels  11,  920  (1997). 

6.  M.  Kotelyanskii,  Trends  Polym.  Sci.  5,  192  (1997). 

7.  R.  Khare,  M.  E.  Paulaitis,  S.  R.  Lustig,  Macromolecules  26,  7203  (1993). 

8.  R.  F.  Rapold,  U.  W.  Suter,  D.  N.  Theodorou,  Macromol.  Theory  Simul.  3, 19-4  (1994). 

9.  R.  J.  Roe,  D.  Rigby,  J.  Phys.  Chem.  89,5280(1988). 

10.  H.  Furuya,  M.  Mondello,  H.  J.  Yang,  R.  J.  Roe,  Macromolecules  27,  5674  (1994). 

11.  R.  J.  Roe,  Computer  Simulation  of  Polymers.  R.  J.  Roe,  Ed.,  Polymer  Science  and 
Engineering,  Prentice-Hall,  New  Jersey  (1991). 

12.  A.  A.  Gusev,  M.  M.  Zender,  U.  W.  Suter,  Macromolecules  27,  615  (1994). 

13.  A.  K.  Rappe.;  W.  A.  Goddard,  III;  Journal  of  Physical  Chemistry  95,  3358  (1991). 

14.  (a)  D.W.  van  Krevelen;  Properties  of  Polymers:  Their  Correlation  with  Chemical  Structure; 
their  Numerical  Estimation  and  Prediction  from  Group  Contributions,  Elsevier  Science 
Publishers,  NY,  pg  536.  (b)  Ibid  p.  575  (1990). 

15.  Belmares,  M.P.,  "Molecular  origins  of  the  thermophysical  properties  of  polymers  and 
modeling  of  polymer  permeation  by  large  molecules",  Caltech  Thesis  (Ph.  D.),  UM 
#9916173,  1998.  The  procedure  has  been  coded  as  a  single  application  using  the  Software 
Developer’s  Kit  (SDK)  in  Cerius2,  Version  4.0,  Accelrys,  Inc.  San  Diego,  CA. 

16.  Originally  part  of  the  suite  of  programs  in  Polygraf,  the  amorphous  builder  is  part  of  the 
Cerius2  software  package,  Cerius2,  Accelrys,  Inc.  San  Diego,  CA. 

17.  S.L.  Mayo,  B.D.  Olafson,  W.A.  Goddard  III  J.  Phys.  Chem.  94,  8897  (1990). 

18.  K.L.  Hoy,  J.  Paint  Technology  42,  79  (1970). 

19.  "Polymer  Handbook",  ed  J.  Brandrup,  E.H.  Immergut,  and  E.  A.  Grulke,  4th  ed.,  J.  Wiley, 
New  York  (1999). 

20.  60th  edition,  CRC  Handbook  of  Chemistry  and  Physics,  ed.,  Robert  C.  Weast  and  Melvin  J. 
Astle,  CRC  Press,  Boca  Raton,  FI,  pgs.  732-735  (1980). 

21 .  Principles  of  Polymer  Systems,  ed.,  Ferdinand  Rodriguez,  2nd  Ed,  McGraw-Hill,  New  York, 
pg  25-26  (1982). 

22.  Physical  Properties  of  Polymers  -  Prediction  and  Control,  A.  Askadskii,  Gordon  and  Breach, 
Amsterdam,  pg.  41-48  (1996). 

23.  2nd  edition,  CRC  Handbook  of  Solubility  Parameters  and  Other  Cohesion  Parameters,  A.  F. 
M.  Barton  Ed.  Boca  Raton,  FI,  pgs.  94-110, 153-159,  CRC  Press  (1991). 

24.  "Graphic  Analysis  of  Resin  Solubilities",  ed.  Jean  P.  Teas  40,  No.  516  (1968). 

25.  Computational  results  obtained  using  software  programs  from  Accelrys,  Inc.  San  Diego,  CA. 
Property/structure  solubility  parameters  calculated  employing  Synthia  program.  Molecular 
dynamics  results  obtained  employing  the  COMPASS  forcefield  with  Amorphous 
Cell/Discover  Programs. 

26.  a.  B.  E.  Eichinger,  D.  R.  Rigby,  and  M.  H.  Muir,  Comp.  Pol.  Sci.  5(4),  147  (1995);  b.  B.  E. 
Eichinger,  D.  R.  Rigby,  and  J.  Stein,  Polymer  43,  5999  (2002);  c.  H.  Sun,  J.  Phys.  Chem.  B. 
102, 7338(1998). 

27.  G.  Walker,  N.  Lewis,  private  communication. 


28.  E.J.  Severin,  B.J.  Doleman,  andN.S.  Lewis,  Analytical  Chemistry  72,  658  (2000). 

29.  M.C.  Lonergan,  E.J.  Severin,  B.J.  Doleman,  S.A.  Beaber  SA,  R.H.  Grubb,  N.S.  Lewis, 
Chemistry  of  Materials  8,  2298  (1996). 

30.  Newman,  S.,  Polymer  Blends,  D.  R.  Paul,  Ed.,  Academic  Press,  Inc.,  Orlando,  vol.  1,  (1978). 

31.  D.  W.  V.  Krevelen,  Properties  of  Polymers :  Their  Correlation  with  Chemical  Structure;  their 
Numerical  Estimation  and  Prediction  from  Group  Contributions,  Elsevier  Science  Publishers, 
New  York,  (1990). 

32.  S.  J.  Plimpton,  R.  Pollock,  M.  Stevens,  Proc  of  the  Eighth  SIAM  Conference  on  Parallel 
Processing  for  Scientific  Computing,  Minneapolis,  MN,  March  (1997). 

33.  S.  J.  Plimpton,  "Fast  Parallel  Algorithms  for  Short-Range  Molecular  Dynamics",  J  Comp 
Phys  117,  1  (1995). 

34.  Severin,  E.J.;  “ Array-based  vapor  sensing  using  conductive  carbon  black-polymer  composite 
thin  film  detectors”,  Caltech  Dissertation,  1999,  ISBN  0-599-42672- 

35.  A.  K.  Rappe,  C.  J.  Casewit,  K.  S.  Colwell,  W.  A.  Goddard  III,  and  W.  M.  Skiff;  J.  Am.  Chem. 
Soc.  114,  10024(1992) 


Table  I.  Comparison  of  calculated  and  experimental  solubility  parameters  for  65  common  solvents,  reactants,  and  monomers  in  order 
of  increasing  experimental  Hildebrand  solubility  parameter. 
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Table  II.  Example  of  output  from  the  CED  molecular  dynamics  method.  Here  the  simulation 
procedure  used  ten  samples  to  estimate  the  condensed  phase  properties  of  2-ethylhexylacrylate. 

Sample  Cohesive  Solubility  Density  End-End  Radius 

Hansen  Solubilities 

Energy  Parameter  distance  Gyration  Elec 

Dispersion  H-Bond 


(cal/cc)  (cal/cc)  1/2  (g/cc) _ (A) _ (A]_ 


(cal/cc)  1/2 

1  -73.42  8.57 

0.80 

7 . 10 

3.19 

3.52 

7.56 

0.00 

2  -81.04  9.00 

0.88 

6.78 

3.19 

4 . 05 

8.31 

0.00 

3  -67.84  8.24 

0.83 

7 . 17 

3.17 

3.20 

7 . 72 

0.00 

4  -77.06  8.78 

0.83 

7 . 63 

3.26 

3.60 

7 . 84 

0.00 

5  -70.63  8.40 

0.84 

7 . 33 

3.17 

3.60 

7 .86 

0.00 

6  -70.46  8.39 

0.82 

6.67 

3.22 

3.18 

7.72 

0.00 

7  -77.91  8.83 

0.87 

7 . 34 

3.13 

3.42 

8.21 

0.00 

8  -81.90  9.05 

0.87 

7 . 17 

3.19 

4 . 00 

8.15 

0.00 

9  -70.78  8.41 

0.85 

6.82 

3.18 

3.81 

7 . 94 

0.00 

10  -65.89  8.12 

0.82 

6.74 

3.19 

3.96 

7 . 87 

0.00 

Density 

Cohesive  Energy  Density 

Average 
0.84 
-73. 69 

Standard  Deviation 
+-  0.03  (g/cc) 

+-  5.52  (cal/cc) 

Solubility  Parameter 

8 . 58 

+-  0.32 

(cal/cc) 

1/2 

Electrostatic  Hansen 

SP 

17 . 55 
3.63 

+-  0.66 
+-  0.32 

(Mpa)  1/2 
(cal/cc) 

1/2 

Dispersion  Hansen 

SP 

7 . 92 

+-  0.24 

(cal/cc) 

1/2 

Table  II  (Example  Output  of  CED 

Simulation) Continued: 

Hydrogen  Bond  Hansen 

SP 

0.00  + 

-  0.00 

(cal/cc) 

1/2 

Non-Bond  EEX 

Unit  Cell  Volume 
End-to-End  Distance 

-76.05 

5820.43 

7 . 07 

+-  5.49  (cal/cc) 

+-174.32  A3 
+-  0.3134  (A) 

Radius  of  Gyration 


3.19  +-  0.0334  (A) 


Table  III.  Average  error  of  Synthia  predictive  methods  compared  with 
CED  predictions  of  Hildebrand  Solubilities. 


Method 

RMS  (a) 

CED-ESP 

0.806 

Synthia-Fedor 

1.388 

Svnthia-van-Krevelen 

1.202 

(a)  Root  mean  square  deviation  (hildebrands). 


Table  IV.  Pearson’s  correlation  coefficients  and  slopes  of  predicted  versus  experimental 
changes  in  resistivity,  AR/R,  for  each  of  seven  polymer  vapor  solvent  detectors. 

Polymer  Sensor _ Slope  Pearson’s  R 


Polycaprolactone 

0.858 

0.925 

Polysulfone 

0.932 

0.962 

PMMA 

0.678 

0.827 

PEVA 

0.888 

0.936 

Polyethylene 

0.870 

0.933 

Polyethyleneoxide 

0.746 

0.874 

Polv(4-hvdroxvstvrene) 

1.018 

0.991 

Table  V.  Estimated  molar  heats  of  vaporization  and  Hansen  solubilities 


A  fi 

*1  vup 

<5  , 

*  2 

. . . 6  3  _J 

Odorants 

Kcal/mol 

Electrostatic  Dispersion 

H-bonding; 

,2-pentanol 

-151.42 

-53.32 

-76.48 

-21.62  1 

3-pentanol 

-142.40 

-47.89 

-76.87 

-17.64 

amylacetate 

-127.31 

-40.19 

-87.13 

0.00 

butylacetate 

-132.03 

-41.75 

-90.28 

0.00  i 

decylacetate 

-104.70 

-21.02 

-83.68 

o.oo  ; 

ethanol 

-257.64 

-146.00 

-51.35 

-60.29  [ 

lethylacetate 

-159.31 

-68.99 

-90.33 

0.00 

\hexylacetate 

-122.55 

-34.83 

-87.72 

0.00 

iso-amylalcohol 

-159.46 

-59.82 

-73.87 

-25.77  | 

\isoamylacetate 

-125.90 

-38.67 

-87.24 

0.00 

isoamylbenzoate 

-119.56 

-23.04 

-96.52 

0.00  : 

isoamylbutyrate 

-111.52 

-25.34 

-86.17 

0.00 

isoamylcaproate 

-104.57 

-20.83 

-83.74 

0.00 

isoamylpropionate 

-113.17 

-30.36 

-82.81 

0.00  [ 

isobutyiacetate 

-130.92 

-45.05 

-85.87 

0.00  ) 

isopropylacetate 

-143.46 

-57.20 

-86.26 

0.00  j 

rn-amylalcohol 

-159.42 

-59.53 

-75.46 

-24.44  | 

n-heptanol 

-130.23 

-37.63 

-76.59 

-16.01  j 

n-hexanol 

-141.38 

-46.42 

-77.97 

-16.99  [ 

n-propanol 

-193.82 

-94.68 

-60.77 

-38.37 

octanol 

-127.59 

-33.80 

-79.91 

-13.88 

octylacetate 

-112.37 

-26.42 

-85.95 

0.00 

propylacetate 

-142.96 

-54.90 

-88.06 

0.00 

n-butanol 

-152.72 

-64.31 

-61.58 

-26.82  ( 

Polymer  Sensor 

PMMA 

-90.51 

-31.19 

-59.32 

0.00  1 

P4HS 

-106.66 

-28.66 

-64.48 

-13.51  j 

PEO 

-168.10 

-68.36 

-95.90 

-3.84 

PE 

-85.45 

-1.00 

-84.46 

0.00 

PEVA 

-85.02 

-10.82 

-74.20 

0.00 

Polysulfone 

-138.74 

-29.76 

-108.98 

0.00 

Caprolactone 

-122.66 

-35.31 

-87.34 

0.00 

Figure  1.  A  polymer  or  solvent  sample  is  put  through  a  series  of  compression  and  expansion 
steps  until  the  proper  density  and  packing  is  obtained.  On  the  left  the  initial  density  is  0.4  pt,, 
40%  of  the  target  density.  After  compression,  second  step,  the  sample  is  over  compressed  to  1 .2 
p0  .  Finally  the  sample  is  allowed  to  relax.  Through  NPT  molecular  dynamics  a  final  prediction 
of  the  density  and  cohesive  energy  of  the  sample  is  obtained.  The  process  is  repeated  for  a  few 
samples  to  gather  statistics. 


Figure  2.  CED  versus  experimental  Hildebrand  solubility  parameters  for  all  molecules  in  Table 
I.  Error  bars  indicate  one  experimental  and  simulation  standard  deviation.  Charge  assignment 
method  is  HF  6-3 1G**  Mulliken  population  charges. 


Figure  3.  CED  versus  experimental  Hildebrand  solubility  parameters  with  quantum  mechanical 
electrostatic  potential  (ESP)  HF  6-3 1G**  assigned  charges.  Error  bars  indicate  one  experimental 
and  simulation  standard  deviation. 


Experimental  AR/R 


Figure  4.  Comparison  between  theory,  equation  (5),  and  experimental  changes  in  resistivity  of 
seven  polymer  sensors  exposed  to  twenty-four  solvents. 


Dr.  Gordon  Shepherd 


Our  subproject  carried  out  several  studies  during  the  five  years  of  the  grant  plus  the  final  one 
year  no-cost  extension  (overall  1998-2004).  The  overall  aim  of  these  studies  was  to  provide 
computational  analysis  of  key  aspects  of  experimental  data  generated  by  our  laboratory  and  by 
our  collaborators  that  would  give  insight  into  the  molecular  and  cellular  basis  of  information 
processing  of  olfactory  sensory  input  in  the  olfactory  pathway. 

Study  1.  Interactions  between  odor  molecules  and  receptor  binding  pocket.  A  major  goal 
of  our  subproject  was  to  gain  insight  into  the  nature  of  the  information  that  is  carried  by  the  odor 
molecules  and  transferred  into  the  neural  domain  by  interactions  with  olfactory  receptor  proteins. 
The  first  study  was  carried  out  by  Michael  Singer,  a  graduate  student  in  our  laboratory.  For  his 
thesis  work,  supported  in  part  by  the  MURI,  tested  by  molecular  modeling  the  ligand  specificity 
of  a  specific  receptor,  17,  that  had  been  shown  in  expression  experiments  to  interact  preferably 
with  octanal,  a  straight-chain  8  carbon  aldehyde.  Using  commercial  software,  he  constructed 
model  of  17  and  carried  out  docking  experiments,  which  showed  that  the  model  also  preferred 
octanal  over  neighboring  aldehydes.  This  analysis  pointed  to  specific  interactions  between  the 
OH  group  of  the  aldehyde  and  the  positive  nitrogen  group  on  a  specific  lysine  residue  in  a 
binding  pocket,  as  well  as  weak  interactions  with  other  residues  forming  the  walls  of  the  pocket. 
This  was  among  the  first  evidence  for  the  specificity  of  the  binding  pocket  in  an  expressed 
olfactory  receptor.  It  was  reported  in  Chemical  Senses  in  2000  (Singer,  2000). 

In  a  parallel  study  we  collaborated  with  the  Goddard  group  in  the  MURI  program  to  carry 
out  a  much  more  extensive  study  of  receptor  binding  properties  in  a  set  of  receptors  that  had 
been  characterized  physiologically  by  the  Buck  laboratory  and  shown  to  be  correlated  with 
specific  receptors  as  identified  by  single  cell  PCR.  This  project  is  described  in  more  detail  in  the 
Goddard  progress  report.  This  three  way  collaboration  resulted  in  a  paper  (Floriano  et  al,  2000) 
which  also  showed  that  the  model  gave  rank  ordering  of  ligand  preference  similar  to  that  seen  in 
the  experiments.  Together,  these  studies  supported  by  MURI  have  provided  some  of  the  first 
and  most  convincing  evidence  for  a  binding  pocket  in  the  7TM  olfactoiy  receptors,  and  for  the 
differential  interactions  between  determinants  of  the  odor  molecules  and  specific  residues  within 
the  binding  pocket.  They  thus  give  direct  insight  into  the  elementary  sensory  units  that  are  the 
basis  for  subsequent  processing  by  circuits  in  the  olfactory  pathway. 

In  subsequent  work  Dr.  Chiquito  Crasto,  with  the  assistance  of  Peter  Bail,  has  continued 
Singer's  work  in  collaborating  with  the  Goddard  and  Buck  laboratories  to  extend  their 
computational  study  to  further  olfactory  receptors  and  their  interactions  with  odor  ligands.  Under 
the  supervision  of  Drs.  Shepherd,  Miller  and  Nadkarni,  he  is  also  responsible  for  the  olfactory 
receptor  database  (ORDB)  and  the  odor  ligand  database  (ODORDB),  funded  by  a  Human  Brain 
Project  grant,  which  support  the  analysis  of  odor  ligand-olfactory  receptor  interactions,  from 
which  our  collaborative  study  has  benefitted. 

It  is  widely  assumed  that  genes  map  directly  into  behavior,  and  therefore  that  the  lower 
numbers  of  olfactory  genes  in  humans  compared  with  rodents  reflect  a  diminished  importance  of 
smell  in  humans.  We  have  argued  that  this  direct  mapping  cannot  be  assumed,  and  that  other 
factors  must  be  taken  into  account,  which  suggest  that  human  smell  is  more  important  than 
generally  realized  (Shepherd,  2004). 

Study  2.  Gene  chip  analysis  of  expressed  olfactory  receptors.  Although  genomic  DNA 
indicated  the  numbers  of  olfactory  receptor  genes  in  human  and  mouse  (350  and  1200 
respectively)  it  left  unknown  how  many  of  these  are  actually  expressed  in  the  olfactory 


epithelium.  To  answer  this  we  have  undertaken  a  collaborative  study  with  the  laboratory  of 
Stuart  Firestein  in  Columbia  and  Minghong  Ma  in  my  laboratory,  who  now  is  at  U  Penn, 
building  the  first  Affymetrix  gene  chip  for  this  purpose.  Thus  far  the  chips  have  been  able  to 
identify  over  800  expressed  genes  (Zhang  et  al,  in  preparation).  Further  studies  will  use  the 
chips  for  analysing  expression  in  different  organs  and  at  different  stages  of  development.  A  new 
database  for  receiving  this  chip  data  is  in  development,  under  the  supervision  of  Drs.  Shepherd, 
Miller  and  Nadkarni. 

Study  3.  Differential  odor  responses  in  populations  of  olfactory  receptor  cells. 

Most  of  our  knowledge  of  the  differential  responses  of  olfactory  receptor  cells  to  different  odors 
has  been  obtained  from  recordings  from  single  cells.  Because  of  the  enormous  number  of 
different  receptors,  a  high  throughput  recording  method  is  needed.  To  answer  this  need  we 
developed  an  isolated  preparation  of  a  segment  ("swatch")  of  the  epithelium,  in  which  we  could 
use  Ca  sensitive  dyes  to  image  the  odor  responses  of  up  to  several  hundred  cells  at  one  time. 

This  enabled  us  to  characterize  for  the  first  time  the  simultaneous  responses  of  many  cells  to  the 
same  odor,  and  to  a  battery  of  different  odors  (Ma  and  Shepherd,  2000).  This  same  approach  has 
been  used  to  analyse  for  the  first  time  the  response  selectivity  of  cells  in  the  septal  organ  of  the 
rat  (Ma  et  al,  2003).  To  extend  these  studies  we  have  collaborated  with  the  laboratory  of  Peter 
Mombaerts  to  use  gene  targeted  mice  with  cells  expressing  the  MOR23  olfactory  receptor  gene 
tagged  with  GFP.  A  report  is  in  preparation  (Grosmaitre  et  al,  2003). 

Study  4,  Differential  odor  responses  in  populations  of  olfactory  glomeruli.  The 
differential  odor  responses  of  receptor  cells  are  projected  onto  the  glomeruli  of  the  olfactory 
bulb,  giving  rise  to  spatial  patterns  of  activity  reflecting  the  differential  responses  of  the 
individual  glomeruli  (also  called  odor  maps  or  odor  images).  Building  on  our  study  of  aldehyde 
responses  of  the  17  receptor  (see  above)  we  have  carried  out  a  study  of  the  odor  images  elicited 
by  the  different  aldehydes.  For  this  we  have  used  high  resolution  fMRI,  which  we  introduced 
several  years  ago  for  olfactory  studies  in  collaboration  with  the  imaging  group  at  Yale  under 
Robert  Shulman  and  Douglas  Rothman.  For  this  study  we  drew  on  our  Human  Brain  Project 
expertise  to  develope  methods  for  rendering  the  3D  image  as  flat  maps  (Liu  et  al,  2004).  The 
results  (Xu  et  al,  2003)  show  that  the  images  are  overlapping  but  different  for  the  different 
carbon  chain  lengths  of  the  aldehydes.  These  fMRI  maps  have  the  advantage  of  being  obtained 
in  the  same  animal  and  of  showing  the  entire  global  extent,  which  is  not  possible  to  obtain  with 
other  methods  (Xu  et  al,  2000).  We  have  pointed  out  the  advantages  of  the  olfactory  bulb  in 
analysing  the  relation  between  activity  and  metabolic  requirements  underlying  fMRI  imaging  in 
the  brain  (Shepherd,  2003). 

Study  5.  Complex  processing  by  neuronal  dendrites  in  the  olfactory  bulb.  The 
information  contained  in  the  odor  maps  is  processed  by  the  microcircuits  of  the  olfactory  bulb. 
These  microcircuits  are  largely  formed  by  the  dendrites  of  mitral/tufted  cells  and  their 
interactions  with  two  types  of  interneuron:  periglomerular  (PG)  cells  and  granule  cells. 
Previously  we  have  shown  that  the  site  of  action  potential  initiation  shifts  dynamically  dependent 
on  levels  of  synaptic  excitation  and  inhibition.  To  test  this  model  further  we  have  carried  out  a 
computational  study,  which  has  shown  a  novel  "ping-pong"  effect  of  the  action  potential 
initiation  site  shifting  from  distal  to  proximal  dendrite  under  different  conditions  of  synaptic 
excitation  in  the  glomerulus  and  synaptic  inhibition  from  the  granule  cells  (Chen  et  al,  2002). 
Further  analysis  of  the  properties  of  the  distal  dendrites  in  the  glomerular  tuft  have  been  carried 
out  using  Ca  imaging  combined  with  two-photon  microscopy  (Zhou  et  al,  in  preparation). 


Mitral  cells  are  also  interconnected  in  their  distal  dendrites  in  the  glomerulus  by  gap 
junctions  which  mediate  electrical  interactions  between  them.  It  has  been  postulated  that  this 
provides  a  mechanism  for  synchronization  of  mitral  cells.  Experimental  evidence  for 
synchronization  has  been  obtained  in  recordings  from  mitral  and  periglomerular  cells  (Zhou  et  al, 
in  preparation;  Jia  et  al,  in  preparation).  To  test  the  gap  junction  hypothesis  we  have  carried  out 
computational  studies,  which  have  shown  the  operating  range  over  which  the  gap  junctions  are 
likely  to  bring  about  synchronization.  The  study  has  also  provided  clear  evidence  of  the  critical 
role  of  the  long  primary  dendrite  in  the  integrative  actions  of  the  mitral  cell,  which  need  to  be 
taken  into  account  in  the  construction  of  neural  networks  of  these  and  similar  cortical  cells.  A 
report  is  in  preparation  (Migliore  and  Shepherd,  in  preparation). 

Study  6.  Neuronal  functional  phenotype  depends  on  dendritic  properties.  The 
mitral/tufted  cells,  like  other  neuron  types,  have  a  distinctive  morphology  and  mode  of  action 
potential  firing.  In  order  to  gain  insight  into  the  relation  between  the  morphological  types  and 
firing  patterns  we  carried  out  a  study  to  classify  them.  The  tentative  result  was  a  flow  diagram 
showing  how  different  types  of  dendritic  morphology  and  different  types  of  firing  patterns  could 
be  related  to  the  different  distributions  of  voltage-gagted  ion  channels  in  the  dendrites  by  a 
hierarchical  scheme  (Migliore  and  Shepherd,  2003).  Further  studies  will  be  needed  to  test  this 
hypothesis. 


Dr.  James  Bower: 


Our  subproject  carried  out  several  studies  during  the  five  years  of  the  grant  plus  the  final  one 
year  no-cost  extension  (overall  1998-2004).  In  general  our  efforts  were  focused  on  two  main 
areas:  physiological  and  modeling  studies  of  the  relationships  between  the  periodic  oscillatory 
behavior  of  the  olfactory  system  and  an  effort  to  understand  how  humans  classify  and  organize 
olfactory  information. 

A.  Is  there  a  natural  order  in  olfactory  stimuli? 

Odor  quality  is  thought  to  depend  primarily  on  the  chemical  structure  of  the  olfactory  stimulus. 

In  form,  this  hypothesis  can  be  traced  back  to  Lucretius,  who  in  50  BCE  suggested  that  agreeable 
odors  were  produced  when  the  olfactory  system  encountered  "smooth"  particles,  while 
disagreeable,  harsh  odors  were  generated  by  "hooked"  particles  [1],  This  idea  was  made  more 
explicit  in  the  early  20th  century  when  Henning  [2]  linked  odor  classes  to  specific  molecular 
features  (as  referenced  in  Schiffman  [3]).  Current  thinking  suggests  that  the  olfactory  system 
functions  as  a  “chemical  classifier”,  where  olfactory  perception  is  organized  around  the  detection 
of  different  structural  classes  of  molecules  by  different  classes  of  receptor  proteins  (e.g.,  [4-7]). 

In  support  of  this  idea,  the  organizational  principles  governing  molecules  in  organic  chemistry 
have  become,  in  practice,  the  favored  hypothetical  organizational  principles  for  olfactory 
perception.  Specifically,  structural  homology  is  implicitly  used  to  explore  odor  “topography”  in 
the  olfactory  system  with  most  studies  based  on  sets  of  chemical  stimuli  representing  small 
changes  in  chemical  structure  between  molecules  with  similar  chemical  properties  (e.g.,  carbon 
chain  length  [8-10]).  Thus  homologous  series  of  alcohols  [8,  11-15],  aldehydes  [8,  11,  16,  17], 
and/or  other  simple  organic  molecules  [18-20]  are  regularly  used  to  examine  the  response 
properties  of  olfactory  receptor  molecules  [20],  olfactory  receptor  neurons  [12],  bulbar  [11, 16, 
19]  and  cortical  [20]  neuronal  responses,  as  well  as  psychophysical  responses  to  these  stimuli  [8, 
13-15]. 

While  probing  the  olfactory  system  with  homologous  series  of  odorants  is  standard  across  a  wide 
range  of  studies,  the  results  obtained  usually  indicate  some  more  complex  pattern  of 
organization.  Thus,  for  example,  studies  of  single  olfactory  receptor  binding  properties  indicate 
that  while  individual  receptors  might  bind  to  two  to  three  consecutive  molecules  in  a  series, 
molecules  taken  from  other  series  also  often  elicit  responses  [20].  More  generally, 
psychophysical  analysis  of  odor  perception  demonstrates  that  similarity  in  molecular  structure 
does  not  also  predict  similarity  in  odor  quality,  an  observation  consistently  observed  for  almost 
half  a  century  [21].  In  fact,  it  is  still  not  possible  to  predict  the  odor  quality  elicited  by  a  stimulus 
based  on  molecular  structure  alone  [3]or  even  to  know  for  certain  that  a  particular  molecule  will 
elicit  an  odor  at  all.  Finally,  the  strongest  argument  that  the  organization  of  odor  perception  is  not 
based  simply  on  structural  homology  is  that  a  change  in  the  concentration  of  the  same 
monomolecular  stimulus  can  produce  dramatic  changes  in  its  elicited  odor  [22],  Thus,  it  is  clear 
that  the  olfactory  system  is  doing  something  more  complex  than  merely  detecting  different 
molecules  based  on  their  chemical  structure. 

In  this  paper  we  report  the  results  of  an  effort  that  initially  disregarded  the  chemical  structure  of 
odorants,  and  instead  sought  evidence  of  order  in  the  relationships  between  odor  perceptions,  and 
subsequently  whether  a  similar  order  might  exist  in  the  organization  of  olfactory  stimuli.  First, 


using  two  large  published  data  sets  of  odor  descriptors  elicited  by  different  monomolecular 
stimuli  [23,  24],  we  have  used  a  cross-entropy  analysis  of  the  co-occurrence  of  pairs  of 
descriptors  associated  with  different  molecules,  to  generate  directed  graphs  representing  the 
relationships  between  these  descriptors.  We  have  found  that  both  datasets,  generated  by 
different  olfactory  specialists,  are  quite  similar.  Second,  we  have  further  analyzed  these  graphs 
with  respect  to  molecules  known  to  elicit  different  odor  perceptions  with  varied  concentration. 
This  analysis  has  revealed  that  all  tested  molecules  of  this  class  fall  into  one  of  three  transitions 
in  the  graph.  Finally,  when  the  molecular  structure  of  the  molecules  whose  descriptors  also  map 
to  these  graph  locations  are  analyzed,  we  discovered  that  the  entire  graph  can  be  divided  into 
regions  whose  descriptors  represent  molecules  containing  sulfur,  nitrogen  or  carbon.  We 
interpret  this  result  to  imply  that  olfactory  perception  reflects  some  understanding  of  different 
biogenetic  pathways  in  the  natural  world,  rather  than  the  kind  of  organization  suggested  by 
structural  chemical  homology.  This  hypothesis  has  important  implications  for  the  interpretation 
of  studies  of  the  structure  and  specificities  of  olfactory  receptors  as  well  as  investigations  of 
structure-function  relationships  in  olfactory  systems  as  a  whole. 
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B.  Relations  between  sniffing  and  oscillations  in  olfactory  cortex 

While  it  has  been  known  for  many  years  that  the  nervous  system  generates  oscillatory  behavior 
in  a  wide  range  of  frequencies  (Adrian  1942,  Freeman  1960,  Bressler  1990,  Steriade  1997),  the 
functional  significance  of  this  dynamical  behavior  has  recently  become  a  growing  focus  for 
many  physiologists  (Engels  and  Singer  2001),  modelers,  and  theorists.  The  oscillatory  behavior 
of  cerebral  cortex  and  its  associated  structures  in  particular  is  now  regarded  by  many  as  being 
directly  linked  to  function  (Engels,  2001).  Current  speculations  on  these  functions  range  from 
the  coding  of  sensory  information  (Singer  and  Gray  1995),  to  a  substrate  for  attention  (Fries 
1997,  Tiitinen  1993,  Tallon-Baudry  1999),  awareness  and  conciousness  (Engels  and  Singer 
2001,  Sauve  1999,  Crich  and  Koch  1990). 

Central  in  many  interpretations  of  the  functional  significance  of  cerebral  cortical  oscillatory 
behavior  is  the  question  of  the  origins  of  this  behavior,  and,  in  particular,  its  relationship  to 
afferent  input.  Modeling  efforts,  including  our  own  (Wilson  and  Bower  1992),  together  with 
experimental  works  (Sanchez-Vives  and  McCormick  1999),  have  suggested  that  the  tendency  of 
cortical  networks  to  oscillate  is  an  intrinsic  property  of  both  its  networks  and  neurons  (Wilson 
and  Bower  1992,  Hasselmo  2001).  Others  have  observed  that  oscillatory  patterns  could  be 
“driven”  by  deeper  brain  structures  like  the  thalamus  (Steriade,  Destexhe)  and  that,  under 
specific  conditions,  they  could  be  related  to  sensory  inputs  (Ferster).  Additional  confusion  in  the 
literature  results  from  the  fact  that  the  cerebral  cortex  generates  many  different  frequencies  of 
behavior,  and  it  is  often  unclear  how  the  oscillations  described  by  one  laboratory,  relate  to  the 
oscillations  described  by  another. 


In  several  studies  funded  by  the  MURI  and  now  complete,  we  have  focused  on  oscillations  in 
the  olfactory  system  occurring  at  frequencies  less  than  1.5  Hz.  These  low  frequency  oscillations 
have  been  recorded  here  under  conditions  of  ketamine/xylazine  anesthesia  and  are  comparable  to 
similar  slow  oscillations  observed  in  neocortical  areas  by  Steriade  and  colleagues  using  the  same 
anesthetic  protocol  (Steriade  1993  a,  b).  Based  on  similarities  between  these  oscillations  and 
those  seen  during  slow  wave  sleep  (SWS),  Steriade  has  proposed  that  the  presence  of  these  slow 
oscillations  reflects  a  behavioral  condition  in  which  the  brain  is  mostly  closed  to  the  external 
environment  and  running  on  its  own  (Steriade  2000,  Timofeev  1996).  This  assertion  is  consistent 
with  electrophysiological  results  demonstrating  slow  wave  oscillations  in  neocortical  slices 
(Sanchez-Vives  2000)  and  slabs  (Timofeev  2000)  which,  of  course,  are  devoid  of  afferent 
sensory  input.  Recently  similar  oscillations  have  also  been  shown  to  occur  in  thalamic  slices 
lacking  afferent  input  (Hughes  2002), 

We  have  used  this  ketamine/xylazine  preparation  to  investigate  the  relationship  between  slow 
oscillations  in  primary  olfactory  (piriform)  cortex  and  the  presence  of  oscillatory  activity  in  its 
afferent  inputs.  Unlike  sensory  neocortex  in  fact  the  overall  temporal  pattern  of  afferent  inputs  to 
the  olfactory  cortex  is  rhythmical  and  controlled  by  the  sniffing  (bressler,  others).  Using 
intracellular  and  extracellular  recording  techniques,  we  demonstrate  here  for  the  first  time  the 
presence  of  ketamine/xylazine-induced  slow  oscillations  in  membrane  potential  and  local  field 
potentials  recorded  in  vivo  from  piriform  cortex.  Further  we  show  that  the  occurrence  of  these 
neuronal  activity  patterns  is  directly  correlated  with  the  natural  breathing  cycle  of  the  rat,  and 
therefore  is  likely  to  be  directly  related  to  periodic  patterns  of  afferent  input  linked  to  respiration. 
In  supporting  to  this  interpretation  we  demonstrate  a  considerable  reduction  in  the  amplitude  and 
regularity  of  oscillations  when  air  bypasses  the  olfactory  system  in  tracheotomized  preparations. 
Pulsing  air  into  the  nostrils  in  this  preparation  restores  the  oscillations,  and  can  directly  entrain 
their  frequency.  We  therefore  conclude  that  there  is  a  direct  relationship  between  slow 
oscillations  in  the  olfactory  cortex  and  in  its  sensory  input.  The  conclusion  that  respiration  is 
reflected  in  the  temporal  behavior  of  piriform  cortex  is  also  consistent  with  studies  conducted  as 
part  of  this  larger  research  project  in  the  olfactory  bulb  of  rats  (Sobel  1993)  and  with  recent 
imaging  results  in  humans  (Sobel  1998). 

We  have  previously  suggested,  based  on  computational  grounds,  that  the  fundamental  neuronal 
architecture  of  cerebral  cortex  may  have  first  evolved  in  the  context  of  the  olfactory  system  and 
then  been  adapted  to  the  use  of  other  sensory  systems  through  the  evolution  of  neocortex 
(Bower,  1995).  We  have  speculated  that  the  presumed  general  usefulness  of  the  type  of 
associative  learning  necessary  for  object  recognition  in  the  olfactory  system  drove  this 
subsequent  development  of  neocortex.  As  a  consequence,  however,  we  suggest  that  neocortex 
was  forced  to  accept  dynamical  patterns  of  behavior  more  clearly  related  to  olfaction  than 
somatosensation,  vision,  audition,  etc.  We  believe  similarities  in  the  neo  and  olfactory  cortical 
mechanisms  for  the  higher  frequency  theta  (7-12  Hz)  and  gamma  (30-50  Hz)  oscillations 
revealed  by  our  realistic  modeling  efforts  (Wilson  and  Bower,  1991;  1992),  support  this  point  of 
view.  With  respect  to  the  slow  frequency  behavior  studied  here,  because  the  olfactory  system  is 
fundamentally  linked  to  the  metabolic  function  of  respiration,  air  continues  to  move  across  the 
olfactory  epithelium  during  the  sleep  state  providing  an  ongoing  periodic  input.  In  neocortex,  in 
contrast,  sleep  is  associated  with  sensory  deafferentation  through  the  action  of  the  thalamus 
whose  decoupled  neurons  in  the  awake  state  become  coupled,  modulating  slow  wave  activity  in 


the  sleep  state  (Bazhenov,  2002).  We  suggest  that  these  slow  wave  thalamic  oscillations 
functionally  correspond  to  the  slow  periodic  activity  generated  by  the  rhythmic  slow  breathing 
associated  with  sleep  in  mammals.  The  essential  point,  and  prediction  is  that  proper  function  of 
cortical  circuitry  during  sleep,  whatever  that  function  is  (Wilson  and  McNaughton,  1993;  Lee 
and  Wilson,  2002),  requires  or  expects  such  slow  periodic  input  because  of  an  evolutionary  link 
between  cortical  circuitry  and  olfaction.  It  remains  to  be  seen  whether  there  is  any  specific 
temporal  relationship  between  rhythmic  breathing  during  sleep,  and  the  timing  of  slow 
oscillations  in  the  thalamus  and  neocortex. 
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A.  Combinatorial  receptor  codes  for  odors 

In  one  part  of  this  project,  we  asked  how  odorant  receptors  (ORs)  in  the  nose  encode  the 
identities  of  a  vast  array  of  environmental  chemicals.  Mice  have  -1000  different  ORs,  each 
expressed  by  a  different  subset  of  olfactory  sensory  neurons  (OSNs)  in  the  nose.  Although  the 
ORs  were  identified  in  1991,  little  was  known  about  the  odorants  recognized  by  individual  ORs 
because  it  was  not  possible  to  functionally  express  ORs  in  heterologous  cells  where  their 
functions  could  be  tested.  To  circumvent  this  problem,  we  developed  an  alternative  strategy  in 
which  we  first  used  calcium  imaging  to  identify  mouse  OSNs  responsive  to  individual  odorants 
and  then  single  cell  reverse  transcriptase  (RT)-PCR  to  identify  the  OR  genes  expressed  in  those 
neurons. 

In  initial  studies,  we  asked  whether  individual  OSNs  actually  do  express  only  one  OR  gene,  as 
predicted.  In  our  experiments  we  used  a  two  step  RT-PCR  protocol.  In  the  first  step  we 
prepared  and  amplified  cDNAs  matching  all  (or  most)  mRNAs  in  the  neuron.  In  the  second  step, 
we  used  OR  primers  to  amplify  cDNAs  from  an  aliquot  of  the  first  PCR  reaction.  Using 
degenerate  primers  matching  conserved  sequence  motifs  in  mammalian  ORs,  we  identified  OR 
genes  expressed  in  >50  single  OSNs  in  these  and  other  experiments.  In  every  case,  we  identified 
only  one  expressed  OR  per  neuron.  We  also  performed  a  number  of  control  experiments  to 
further  investigate  not  only  the  number  of  OR  genes  expressed  per  OSN,  but  also  the  ability  of 
our  single  cell  RT-PCR  method  to  accurately  identify  the  OR  genes  expressed  in  individual 
OSNs.  The  following  results  provided  strong  evidence  that  each  OSN  expresses  only  one  OR 
gene: 

1 .  OR  cDNAs  obtained  from  the  same  OSN  with  different  primers  were  identical  in  overlapping 
regions. 

2.  All  60  cloned  OR  cDNAs  derived  from  the  same  neuron  hybridized  to  the  same  OR  probe  in 
each  case  tested. 

3.  When  total  cDNAs  amplified  from  pairs  of  cells  were  mixed  and  then  subjected  to  PCR  with 
OR  primers,  the  OR  cDNAs  obtained  from  each  alone  were  present,  indicating  that  more  than 
one  OR  cDNA  could  be  detected  if  present. 

4.  When  the  first  strand  cDNAs  made  from  individual  neurons  were  split  into  3  tubes  and  then 
subjected  to  the  two  step  PCR  protocol,  the  OR  cDNAs  amplified  from  the  3  tubes  were 
identical,  excluding  the  possibility  that  the  ORs  detected  with  our  method  derive  from 
contaminating  genomic  DNA  (since  there  are  only  two  alleles  of  each  gene). 

For  test  odorants,  we  used  a  series  of  aliphatic  odorants  with  straight  carbon  chains  ranging  in 
length  from  4  to  9  carbon  atoms.  We  used  four  different  classes  of  aliphatic  odorants  with  the 
same  carbon  chains,  but  different  functional  groups:  carboxylic  acids,  aliphatic  alcohols, 
bromocarboxylic  acids,  and  dicarboxylic  acids.  Only  a  small  percentage  of  neurons  (0-7%  of 
>600  neurons  tested)  responded  to  individual  test  odorants. 


We  identified  ORs  expressed  in  14  neurons  responsive  to  the  test  odorants.  Since  each  neuron 
expressed  only  one  OR  gene,  the  response  profile  of  the  neuron  indicated  the  recognition  profile 
of  the  OR  it  expressed.  Consistent  with  previous  functional  studies  of  OSNs,  we  found  that 
individual  ORs  can  recognize  multiple  odorants.  Our  results  also  showed  that  a  single  odorant 
can  be  recognized  by  multiple  ORs.  Importantly,  however,  different  odorants  were  recognized 
by  different  combinations  of  ORs.  Thus  the  identities  of  different  odorants  are  encoded  by 
different  combinations  of  receptors.  This  combinatorial  use  of  the  OR  family  should  allow,  in 
principle,  for  an  ability  to  detect  and  discriminate  an  enormous  variety  of  odorants  whose 
number  vastly  exceeds  that  of  ORs.  Even  if  each  odorant  were  encoded  by  only  three  receptors, 
this  combinatorial  coding  scheme  would  generate  almost  one  billion  different  odor  codes. 

Our  studies  also  showed  that  there  is  tremendous  diversity  in  the  recognition  properties  of 
different  ORs.  Both  the  carbon  chain  length  of  the  odorants  and  their  functional  groups  were 
important  for  recognition  for  most  of  the  ORs  we  identified,  but  the  combinations  of  carbon 
chain  lengths  and  functional  groups  that  could  be  recognized,  or  tolerated,  by  individual 
receptors  varied  tremendously.  In  addition,  we  found  that  a  single  odorant  can  be  recognized  by 
both  highly  related  and  divergent  ORs. 

In  addition,  these  studies  showed  that  even  a  slight  change  in  the  structure  of  an  odorant,  or  a 
change  in  its  concentration,  can  change  its  "receptor  code",  the  combination  of  ORs  that 
recognize  the  odorant.  This  suggests  that  changes  in  receptor  code  may  underlie  perceptual 
alterations  in  humans  that  can  result  from  changes  in  odorant  structure  or  concentration. 

In  summary,  our  studies  indicate  that  the  olfactory  system  uses  a  combinatorial  receptor  coding 
scheme  to  distinguish  odors.  In  this  scheme,  the  identities  of  different  odorants  are  specified  by 
different  combinations  of  receptors,  but  each  receptor  can  be  used  as  one  component  of  the  codes 
for  many  odors.  This  scheme  explains  how  1000  different  receptors  can  allow  for  the 
discrimination  of  thousands  or  tens  of  thousands  of  different  odors  in  the  external  environment. 

B.  Detection  of  odorants  and  pheromones  by  the  vomeronasal  organ 

In  addition  to  odorants,  the  olfactory  system  detects  pheromones,  chemicals  released  from  animals  that 
stimulate  hormonal  changes  or  specific  behaviors,  such  as  mating  or  aggression,  in  members  of  the  same 
species.  Many  mammals  have  a  second  olfactory  sense  organ,  the  vomeronasal  organ  (VNO),  which  is 
thought  to  be  specialized  to  detect  pheromones.  In  previous  studies,  we  and  others  found  that  the  VNO 
and  nose  use  different  mechanisms  to  detect  sensory  ligands.  Two  distinct  families  of  receptors  are 
expressed  in  the  VNO,  the  V1R  and  V2R  families,  and  the  different  families  are  expressed  by  different 
subsets  of  VNO  neurons.  As  in  the  nose,  however,  it  appears  that  each  neuron  expresses  only  one 
receptor  gene.| 

To  determine  how  the  VNO  system  encodes  sensory  information,  we  used  calcium  imaging  to 
analyze  responses  of  single  mouse  VNO  neurons  to  pheromones.  We  also  examined  responses 
of  these  neurons  to  odorants. 


We  found  that,  similar  to  OSNs,  VNO  neurons  exhibit  transient  increases  in  intracellular  calcium 
when  they  are  exposed  to  140mM  KC1.  These  increases  are  likely  to  be  due  to  influx  of  calcium 
through  voltage-sensitive  calcium  channels  that  open  upon  membrane  depolarization  by  KC1. 

In  initial  studies,  we  tested  VNO  neurons  for  responsiveness  to  six  different  mouse  pheromones: 
dehydro-exo-brevicomin,  2,5-dimethylpyrazine,  E-beta- farnesene,  2-heptanone,  lactol  (6- 
hydroxy-6-methyl-3-heptanone),  and  2-sec  butyl-4, 5-dihydrothiazole.  The  pheromones  were 
tested  at  lOOuM  except  for  dehydro-exo-brevicomin  (7uM),  2,5-dimethylpyrazine  (luM)  and  2- 
sec  butyl-4, 5-dihydrothiazole  (7uM),  because  of  limited  supplies  when  the  experiments  were 
done.  In  later  studies  supported  by  the  NIH,  we  tested  additional  neurons  with  these  pheromones 
as  well  as  with  odorants.  Interestingly,  individual  pheromones  elicited  responses  in  only  about 
0.3-0. 7%  of  VNO  neurons.  Given  that  there  are  about  260  different  types  of  VNO  receptors,  this 
raises  the  possibility  that  some  pheromones  may  be  detected  by  only  a  single  receptor  type  rather 
than  by  a  combination  of  receptors  as  in  the  nose.  Surprisingly,  we  found  that  VNO  neurons  can 
also  detect  some  odorants,  but  not  all.  Given  the  neural  pathway  followed  by  signals  generated 
in  the  VNO,  this  raises  the  possibility  that  some  odorants  may  also  be  capable  of  stimulating 
hormonal  changes  or  instinctive  behaviors  in  mice. 

C.  The  human  OR  gene  repertoire 

In  the  later  part  of  this  project,  we  characterized  the  human  OR  family.  To  identify  human  OR  genes, 
we  searched  the  NCBI  finished  and  unfinished  human  genome  sequence  databases.  In  these  searches, 
we  used  as  queries  either  short  amino  acid  sequence  motifs  conserved  among  mammalian  ORs  or  a 
diverse  set  of  mouse  ORs.  With  these  queries,  we  searched  for  related  sequences  in  all  6  possible 
translated  reading  frames  of  the  human  genome  sequence.  We  then  isolated  the  DNA  sequence 
surrounding  each  match  and  translated  it  to  determine  the  encoded  protein.  Classification  as  an  OR  gene 
was  based  on  visual  inspection  and,  in  some  cases,  by  searching  the  NCBI  protein  database  with  the 
translated  gene. 

In  these  studies,  we  identified  339  intact  OR  genes  and  297  OR  pseudo  genes  in  the  human  genome. 
Since  the  genome  sequence  was  -93%  complete  at  the  time  we  did  this  work,  this  is  likely  to  represent 
the  full,  or  nearly  full,  repertoire  of  human  ORs. 

We  next  determined  the  chromosomal  location  of  each  human  OR  gene.  We  did  this  by  examining  the 
chromosomal  assignment  of  the  clone  in  which  each  gene  was  found.  Remarkably,  we  identified  51 
different  OR  gene  loci.  We  found  OR  gene  loci  on  every  human  chromosome  except  chromosome  8 
and  the  Y  chromosome,  a  total  of  21  different  chromosomes. 

OR  genes  can  be  divided  into  subfamilies  on  the  basis  of  sequence  relationships.  Previous  studies 
support  the  idea  that  closely  related  OR  genes  that  belong  to  the  same  subfamily  may  encode  receptors 
that  recognize  structurally-related  odorants.  Using  a  cutoff  of  60%  nucleotide  sequence  identity  to 
define  members  of  the  same  subfamily,  the  intact  human  ORs  assorted  into  172  different  subfamilies. 
We  found  that  members  of  the  same  subfamily  are  usually  found  at  the  same  chromosomal  locus.  This 
raises  the  possibility  that  different  OR  gene  loci  might  have  different  functions  in  odorant  recognition. 


It  has  been  proposed  that  humans  may  be  relatively  defective  in  olfactory  perception  compared 
to  mice  and  rats.  This  idea  was  based  on  preliminary  indications  that  mice  may  have  more  OR 
genes  than  humans.  In  other  studies,  which  were  supported  by  the  NIH,  we  analyzed  the  mouse 
OR  gene  repertoire.  Comparison  of  human  and  mouse  OR  families  showed  that  mice  have  about 
2.7  times  as  many  ORs  as  humans  and  241  OR  subfamilies  compared  to  the  172  found  in  human. 
When  we  performed  a  detailed  comparison  of  human  and  mouse  ORs,  we  found  that  most 
human  OR  subfamilies  are  found  in  mouse.  Interestingly,  however,  the  mouse  subfamilies 
almost  invariably  have  more  member  ORs.  This  suggests  that  humans  and  mice  are  likely  to 
detect  many  of  the  same  odorant  structural  features,  but  that  mice  may  have  a  superior  ability  to 
discriminate  odorants  with  highly  related  structures. 
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The  process  whereby  an  odorant  (physical  entity)  becomes  a  smell  (percept)  consists  of  several 
stages.  Initially,  the  odorant  must  be  sniffed  into  the  nasal  passage.  This  crucial  initial  stage  of 
olfactory  processing  is  often  overlooked.  Following  this  sampling  phase,  the  odorants  are 
transduced  by  receptors  ligning  the  olfactory  epithelium.  The  transduced  information  is  then 
processed  first  at  the  level  of  olfactory  bulb,  and  later  olfactory  cortex.  This  MURI  investigated 
olfactory  processing  from  Detection  to  Classification,  and  our  lab  contributed  to  the  study  of 
both  processes: 

Odorant  detection  --  Sampling 

Sniffs  are  modulated  in  response  to  odor  content.  Higher  concentrations  of  odor  induce  lesser- 
volume  sniffs.  This  phenomenon  implicates  a  neural  feedback  mechanism  that  measures  sensory 
input  (odor  concentration)  and  modulates  motor  output  (sniffing)  accordingly1,2.  We  are 
interested  in  this  modulatory  mechanism  both  in  order  to  better  understand  the  olfactory  process, 
and  to  explore  the  implementation  of  such  a  mechanism  in  artificial  olfactory  devices  (electronic 
noses).  We  therefore  first  set  out  to  ask  what  are  the  temporal  parameters  of  this  mechanism.  A 
stainless-steel  computer-controlled  olfactometer,  equipped  with  mass  flow  controllers, 
temperature  and  humidity  control,  and  on-line  photo-ionization  detection,  was  coupled  to  a 
highly  sensitive  pneumatotachograph  that  measured  nasal  flow.  The  olfactometer  was  used  to 
generate  four  ascending  concentrations  of  the  odorants  propionic  acid  and  phenethyl  alcohol. 

We  found  that  sniff  volume  was  inversely  related  to  odor  concentration  (p  <  0.0001). 
Furthermore,  sniffs  were  uniform  and  concentration-independent  for  the  initial  150  msec,  but 
acquired  a  concentration-dependent  flowrate  as  early  as  1 60  msec  following  sniff  onset  for 
propionic  acid  (p  <  .05),  and  260  msec  for  phenethyl  alcohol  (p  <  .05).  We  published  these 
results  in  the  MURI-supported  manuscript  by  Johnson  et  al3. 

These  findings  are  important,  firstly  because  they  merit  reframing  our  view  of  the  temporal 
aspects  of  olfaction.  Olfaction  is  traditionally  considered  a  "slow"  sensory  process.  Our  results 
suggest  otherwise,  and  point  to  rapid  processing  in  olfaction.  These  results  also  raise  important 
questions  regarding  olfactory  sensory  transduction.  The  current  view  is  that  odorant  transduction 
takes  around  150  ms4,5.  However,  the  sniff-modulation  that  we  see  within  160  ms,  renders  the 
150  ms  value  unlikely  (it  is  unreasonable  that  the  modulation  would  be  achieved  within  10  ms). 
Our  finding  of  rapid  processing  has  now  also  been  replicated  in  rats6,  suggesting  that  the 
temporal  dynamics  of  transduction  should  be  revisited  and  revised. 

Furthermore,  the  rapid  olfactomotor  function  we  have  characterized  suggests  that  olfactomotor 
sniff  feedback  control  is  subcortical  and  may  rely  on  neural  mechanisms  similar  to  those  that 
modulate  eye  movements  to  accommodate  vision  and  ear  movements  to  accommodate  audition. 
We  have  now  embarked  on  localizing  this  mechanism  in  the  brain,  and  have  some  evidence 
suggesting  that  it  may  be  cerebellar. 


Figure  1.  First  second  of  the  mean  sniff  to  low  and  high  concentrations  of  a)  propionic  acid,  b) 
phenethyl  alcohol.  The  blue  line  is  the  mean  sniff  to  low  concentration,  and  the  red  line  is  the 
mean  sniff  to  high  concentration.  Bars  are  standard  error.  The  p-value  from  the  associated 
paired  t-test  is  shown  in  yellow.  Sniffs  of  propionic  acid  are  significantly  concentration 
dependent  by  160  msec,  and  sniffs  of  PEA  are  significantly  concentration  dependent  by  260 
msec. 


Odorant  detection  --  Specific  anosmia 

The  olfactory  system  is  a  model  for  neural  plasticity  in  that  its  sensory  neurons  and 
connectivity  are  continuously  renewed  throughout  life8.  An  example  of  plasticity  in  the  adult 
human  olfactory  system  is  learned  detection  of  androstenone  (5  alpha-androst-16-en-3-one). 
Approximately  30%  of  the  adult  human  population  does  not  perceive  an  odor  when  sniffing  the 
steroidal  compound  androstenone,  but  such  sensitivity  can  be  induced  by  repeated  exposure7.  To 
understand  the  mechanisms  of  plasticity  underlying  this  acquired  capability,  one  must  first 
determine  its  location  in  the  olfactory  system.  Is  plasticity  occurring  peripherally  at  the  olfactory 
receptors  in  the  nose,  or  centrally  in  the  brain?  To  address  this  question  we  systematically 
exposed  only  one  nostril  of  androstenone  non-detectors  to  androstenone  and  then  tested  the 
unexposed,  naive,  nostril.  Considering  that  the  two  olfactory  pathways  are  not  neurally 
connected  at  the  peripheral  level,  if  the  naive  nostril  had  learned  to  detect  androstenone,  this 
would  suggest  that  plasticity  occurred  centrally.  By  contrast,  if  the  naive  nostril  remained  unable 


to  detect  androstenone  while  the  exposed  nostril  learned,  this  would  suggest  that  plasticity 
occurred  peripherally. 

We  screened  42  subjects  for  androstenone  detection  using  a  four-trial,  three-alternative- 
forced-choice  paradigm  with  criteria  set  at  two  or  fewer  correct  trials.  This  screen  yielded  12 
non-detectors,  a  prevalence  of  non-detection  (29%)  equivalent  to  that  in  previous  reports.  These 
non-detectors  returned  to  the  lab  for  ten-minute  sessions  daily,  for  21  days,  where  one  nostril 
was  continuously  exposed  to  androstenone.  The  other  nostril  was  blocked  by  insertion  of  an 
inflatable  plug.  Five  liters  per  minute  of  heated,  humidified,  bottled  air  were  injected  through  the 
plug  to  prevent  any  androstenone  from  entering  the  occluded  nostril  by  reverse  flow  (retronasal 
olfaction).  At  the  beginning  of  the  experiment  this  group  was  at  25%  accuracy  (SD=1 1%), 
which  is  not  significantly  different  from  33%  chance  (Binomial  from  chance:  p  =  .9).  Twenty- 
one  days  of  exposure  doubled  androstenone  detection  accuracy  in  the  exposed  (change  from 
baseline:  t(l  1)  =  3.3,  p  <  .007.  Binomial  from  chance:  p  <  .002),  and  unexposed  (change  from 
baseline:  t(l  1)  =  2.3,  p<  .04.  Binomial  from  chance:  p  <  .02)  nostrils  respectively  (Fig.  2). 

There  was  no  significant  difference  in  the  extent  of  improvement  between  the  exposed  and 
unexposed  nostrils  (t(l  1)  =  .4,  p  =  .7).  This  finding  indicates  that  the  site  of  plasticity  underlying 
learned  detection  of  androstenone  communicates  with  both  nostrils  via  a  central  component.  We 
published  these  results  in  the  MURI-supported  manuscript  by  Mainland  et  al9 
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Figure  2.  Detection  accuracy  at  pre-exposure  and  in  the  exposed  and  unexposed  nostrils 
following  21  days  of  exposure. 

We  took  advantage  of  the  paired  anatomy  of  the  olfactory  system  to  find  that  the  plasticity 
underlying  the  emergence  of  androstenone  detection  has  a  central  component.  Such  a 
component  may  be  likened  to  pattern  recognition,  Despite  the  equivalent  retinal  image,  a  CT 
scan  may  appear  all  gray  to  the  layperson  while  revealing  intricate  anatomy  to  the  trained 
neurologist.  Similarly,  the  nose  may  have  been  delivering  the  same  signal  to  the  brain  at  the 
beginning  and  end  of  this  experiment,  but  through  repeated  exposure  the  brain  may  have  learned 
to  make  sense  of  a  previously  senseless  signal.  Such  learning  may  have  occurred  at  the  olfactory 


bulbs  or  in  primary  olfactory  cortex,  a  substrate  that  shares  information  from  both  nostrils10  and 
is  optimized  for  olfactory  learning11. 

Odorant  classification 

Two  aspects  of  odor  that  are  encoded  by  the  olfactory  system  are  intensity  (high/low)  and 
valence  (pleasant/unpleasant).  These  aspects  are  tightly  linked.  Most  pleasant  odorants  become 
more  pleasant  as  their  intensity  increases,  and  most  unpleasant  odorants  become  more  unpleasant 
as  their  intensity  increases12.  Some  odorants,  such  as  indol  and  skatol,  dramatically  change  their 
percept  as  intensity  is  varied.  In  other  words,  intensity  is  an  axis  of  odorant  classification.  Here 
we  used  fMRI  to  ask  if  these  tightly  linked  dimensions  of  intensity  and  valence  are  encoded  by 
dissociable  neural  mechanisms.  To  dissociate  the  dimensions  of  intensity  and  valence,  we 
devised  four  stimuli  using  two  odorants.  High  and  low  concentration  versions  of  the  odorants 
citral  and  valeric  acid  were  prepared.  Citral  smells  lemony,  fruity,  and  fragrant,  and  is  perceived 
by  most  as  pleasant13.  Valeric  acid  smells  sweaty,  rancid,  and  sickening  and  is  perceived  by 
most  as  unpleasant13.  The  concentrations  used  were  selected  through  a  psychophysical  prestudy 
where  the  high  concentration  citral  was  rated  similar  in  intensity  to  the  high  concentration  valeric 
acid,  and  the  low  concentration  citral  was  rated  similar  in  intensity  to  the  low  concentration 
valeric  acid.  Further,  an  intensity  range  was  selected  such  that  odor  valence  could  be 
manipulated  with  relative  independence  from  intensity.  Thus,  the  above  design  provided  four 
conditions  [1.  intense-pleasant  2.  intense-unpleasant  3.  unintense-pleasant  4.  unintense- 
unpleasant]  that  allowed  for  a  dissociation  of  intensity  from  valence  (Fig.  3). 
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Figure  3:  Ideal  and  observed  odor  space. 

The  abscissa  represents  odor  intensity.  The  ordinate  represents  odor  valence,  with  increasing 
pleasantness  represented  as  greater  magnitudes.  ( a)  Odor  selections  were  guided  by  an  attempt 
to  construct  an  affective  space  whereby  odor  intensity  and  valence  were  manipulated 
independently.  High  and  low  concentrations  of  the  pleasant-smelling  citral  and  unpleasant¬ 
smelling  valeric  acid  were  psychophysically  pre-selected  to  represent  the  four  critical  design 


quadrants.  (  b)  Subjective  intensity  and  valence  estimates  obtained  after  scanning.  Horizontal 
error  bars  reflect  standard  error  of  the  mean  (s.e.m.)  for  intensity  ratings  from  1  (low  intensity)  to 
9  (high  intensity)  and  vertical  reflects  s.e.m.  for  valence  from  1  (very  unpleasant)  to  9  (very 
pleasant).  Consistent  with  the  manipulation  of  odor  valence,  valeric  acid  was  rated  as 
significantly  more  unpleasant  than  citral  ( FI, 15  =  25.30,  P<  0.0001).  Consistent  with  the 
manipulation  of  odor  intensity,  high-concentration  odors  were  rated  as  significantly  more  intense 
than  low-concentration  odors  ( FI, 15  =  78.45,  P<  0.0001).  Independent  manipulation  of  intensity 
and  valence  was  achieved  for  three  of  the  four  design  quadrants;  this  was  sufficient  to  separate 
neural  responses  to  intensity  and  valence.  Low  concentrations  of  citral  and  valeric  acid  were 
equated  for  intensity  ( FI,  15  =  0.10,  P>  0.75),  but  differed  in  valence  ( FI,  15  =  6.46,  P<  0.02); 
high  and  low  concentrations  of  valeric  differed  significantly  in  intensity  ( FI,  15  =  20.31,  P< 
0.0004),  but  not  in  valence  ( FI, 15  =  0.08,  P>  0.78).  High  concentrations  of  citral  were  rated  as 
more  intense  ( FI,  15  =  10.25,  P<  0.007)  than  high  concentrations  of  valeric,  and  more  pleasant  ( 
FI, 15  =  7.68,  P<  0.01)  than  low  citral.  Consistent  with  their  relative  independence,  there  was  no 
significant  association  between  individuals'  estimations  of  odor  intensity  and  valence  for  both 
odorants  (valeric,  r=  -0.19,  P>  0.27;  citral,  r=  0.29,  P>  0.11). 

We  delivered  the  four  stimuli  interspersed  with  a  clean  air  stimulus  in  a  random  ordered  event- 
related  design  to  16  participants  while  measuring  brain  activation  with  a  3  Tesla  fMRI  scanner. 
We  found  that  amygdala  activation  was  associated  with  intensity  and  not  valence  of  odors. 
Conversely,  distinct  regions  of  orbitofrontal  cortex  were  associated  with  valence  independent  of 
intensity  (Figure  4).  We  published  these  results  in  the  MURI-supported  manuscript  by  Anderson 
et  al14. 

Consistent  with  the  sensory  nature  of  odor  intensity,  in  the  olfactory  system  initial 
intensity  coding  occurs  at  the  levels  of  the  olfactory  epithelium  and  bulb15,  and  has  been 
implicated  at  primary  olfactory  cortex1617,  which  includes  the  amygdala18.  In  contrast,  our 
findings  suggest  that  higher-order  hedonic  differentiation  occurs  in  secondary  olfactory  regions, 
which  include  the  orbitofrontal  cortices18.  Thus,  the  present  results  suggest  a  major  division  of 
the  neural  representation  of  affective  space  into  more  primary  intensity  and  higher-order  hedonic 
components.  In  contrast  with  relatively  preserved  hedonic  responses  in  patients  with  amygdala 
lesions19,  it  has  been  long  known  from  stroke  patients  that  the  prefrontal  cortices  make  important 
contributions  to  the  experience  of  hedonic  tone20’21,  with  more  recent  studies  highlighting 
specifically  the  critical  role  of  the  orbitofrontal  cortices21,22.  Such  critical  prefrontal 
contributions  to  hedonicity  dovetail  nicely  with  appraisal  theories  of  human  emotion  that 
emphasize  how  affective  responses  are  not  simply  a  reflection  of  the  intrinsic  quality  (positive 
vs.  negative)  of  a  stimulus,  but  result  from  interactions  among  the  person,  the  situational  context, 
and  the  stimulus23.  Such  flexibility  of  hedonic  response  is  a  hallmark  of  orbitofrontal 
representations,  which  are  modulated  by  changes  in  affective  relevance,  as  in  the  hungry  versus 
sated  state  of  the  perceiver24.  Thus,  unlike  the  evolutionarily  conserved  functions  of  the 
amygdala,  it  appears  that  the  malleability  of  human  hedonic  experience  is  characteristic  of  the 
more  flexible,  integrative,  and  evolved  functions  of  the  prefrontal  cortices. 


in  ten  sity  Valence  (-/+•) 


Figure  4:  Functional  regions  of  interest  (fROI)  defined  by  their  correlation  with  individual 
differences  either  in  the  evaluation  of  intensity,  pleasantness  or  unpleasantness  of  the  four 
odor  conditions  and  clean  air. 


Scatter  plots  depict  the  degree  of  association  between  individuals'  fROI  signal  and  valence  and 
intensity  evaluations  (depicted  as  circles  in  standardized  units).  The  ordinate  represents  fMRI 
signal,  and  the  abscissa  represents  either  the  evaluation  of  odor  intensity  or  valence  for  each 
stimulus  condition  for  each  participant.  (  a)  A  region  in  the  right  hemisphere  extending  from  the 
dorsal  amygdala  into  the  piriform  cortex  was  correlated  with  the  evaluation  of  intensity  but  not 
valence.  (  b)  A  bilateral  subcallosal  gyrus  activation  extending  into  the  posteromedial 
orbitofrontal  cortex  was  correlated  with  the  evaluation  of  pleasantness  but  not  intensity.  (  c)  Left 
anterior  lateral  and  right  anterior  medial  orbital  cortical  activations  were  correlated  with  the 
evaluation  of  unpleasantness  but  not  intensity. 
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